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PREFACE 


The main objective of this study is to develop a numerical model 
for simulating the movement of well graded sediment through a stream 
network. In Part 2 of this report the theoretical background of the 
model is described and the governing equations are formulated. The 
numerical solution of these equations is presented in Part 3., Hydraulic 
routing can be performed using any acceptable algorithm supplied by the 
user because the water movement is assumed to be uncoupled from the 
sediment process. The model can be used in conjunction with any 
suitable sediment yield model to supply the water and sediment runoff 
from lateral areas. The application and results of the model are 
discussed in Part 4. The computer program for the numerical model of 
the East Fork River system is described and listed in the third 
addendum. 
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INTRODUCTION 


This report describes a one-dimensional numerical model designed to 
simulate sediment transport in natural channels. The physical processes 
associated with the sediment movement are reproduced using a variety of 
algorithms. These algorithms incorporate sets of equations that operate 
on input data in a predetermined sequence to generate output data repro¬ 
ducing the actual physical process. A satisfactory model must include 
the more relevant aspects of that process. 

Sediment moves driven by hydrodynamic forces exerted by the flow of 
water which in many instances is highly time dependent. The sediment 
transport model must, therefore, account for unsteadiness in sediment 
movement. 

The dependence of sediment motion on flow conditions makes it also 
dependent on the longitudinal variations the flow experiences as a 
result of stream boundary irregularities. These variations are 
reflected in the spatial variability of the sediment load distribution. 
Depending on particle size, some particles may be carried primarily in 
suspension, while others move entirely as bed load. In addition, 
depending on flow conditions, particles moving in suspension at one 
place may be moving as bed load farther downstream. ^ Whether a 
particular size fraction moves primarily as suspended load or bed load 
determines to what extent that fraction of the sediment load will lag 
behind the flood wave and therefore determines what the magnitude of 
longitudinal sorting will be. This means that the transport model must 
reflect the dependence of the sediment load lag on the material 
properties of the sediment as well as on hydraulic conditions. Whenever 
the bed material consists of a mixture, its transport involves the 
motion of a multitude of particles of diverse sizes. Some particles may 
deposit on the streambed while others are scoured away, resulting in a 
size composition of the material in transport different from that of the 
bed. A realistic model must account for the interchange between the bed 
surface material and the moving sediment load, and should simulate the 
residual transport capacity of the stream. The latter is a measure of 
the ability of the flow to further entrain material of a given size 
fraction in the presence of all the fractions already in motion. 


During the above particle interchange, the bed material particle 
size composition changes continuously and, in the process, the bed may 
experience a net amount of aggradation or degradation. For certain flow 
conditions the bed degradation in a reach may be limited by the 
formation of an armoring layer, over which sediment may move either in 
suspension or as intermittent bed-load waves. The armoring layer may be 
destroyed during high flows and reformed at subsequent low flows. A 
model must therefore be capable of tracking the streambed profile 
evolution and the changes in bed material size distribution. 

The proposed model is designed to meet the above criteria. It can 
simulate the unsteady transport of sediment mixtures through a network 
of nonbifurcating channel reaches, and it can be used in tandem with a 
suitable sediment yield model, like the one described by Borah et al. 
(1981), which simulates the supply of water and sediment runoff from 
adjacent upland areas. At present the model is restricted to 
consideration of noncohesive bed materials only. It is also restricted 
to down channel streamflow, and cannot consider the effect of transverse 
currents. 

Parts 2 and 3 of this report provide a detailed description of the 
model. Part 4 discusses the validation of the model on sets of 
laboratory and field data. Coding details are given in Addendum 3. 
Some of this material has been presented in an earlier report (Borah, 
1979). 
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2 MODEL FORMULATION 
2.1 EQUATIONS OF MOTION 

The present model treats the time-dependent problem of one¬ 
dimensional sediment routing in alluvial channels. The hydraulic 
functions driving the sediment movement are the local flow discharge, 
stage, and cross-sectional flow area. Therefore, an unsteady flow 
algorithm is required to compute the time and space distribution of 
those flow parameters in the channel, given information concerning chan¬ 
nel geometry, history of inflowing water discharge, and/or downstream 
stage. 

One-dimensional hydraulic routing algorithms are usually based on 
the shallow-water equations of momentum and conservation of mass for 
sediment-laden water. In natural streams load concentrations of up to 
50,000 ppm will not change the density of the mixture by more than .3%. 
Thus, density variations may be ignored. Furthermore, within that range 
of concentrations, the surface waves propagate with a velocity that is 
practically unaffected by the erodibility of the bed (Gradowczyk, 1968). 
The governing equations for water movement can thus be written 
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where A is the flow cross-sectional area, g is the acceleration of 

gravity, Q is the flow discharge, q is the lateral water inflow per 

w 

unit length of channel, S q is the bed slope, is the friction slope, t 
is t i.me, V is the mean flow velocity, is the velocity component of 
lateral inflow in the main flow direction, x is the horizontal distance 
along the channel, y is the flow depth, P is the momentum coefficient, 
and p is the water density (Fig. 1). 

Two approximations to Eqs. 1 and 2 have found wide application in 
unsteady flow routing. They are the diffusion and kinematic wave ap¬ 
proximations (Ponce, Li, and Simons, 1978). The first assumes that the 
inertia terms are negligible, while the second assumes that both the 
inertia and pressure terms can be neglected. Both approximations have 
neen shown to be satisfactory in a variety of cases. Solutions of the 
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complete shal low-water equations and their v.i; n-n.: ippri.x i mat i oris have 
been extensively discussed in the literature (Mi tier and Yevjevich, 
1975) and they will not be considered in this payer Anv acceptable 
hydraulic algorithm will suffice since the water movement is assumed to 
be uncoupled from the sediment processes. 

The third equation i t motion is given by the conservation of mass 
equation for sediment. As shown in Addendum i , this equation can be 
expressed in the form 
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where c is the average volume concentration. Q i s the sediment, volume 
flux, A is the bed poros«* y, W is the active width of the bed (i.e., 
that portion of the fed width in which erosion or deposition takes 
place), z is the local bed elevation, r is the net sediment volume flux 
through the suspended-bed load interface, and q is the lateral inflow 
of sediment per unit channel length. 

In Eq. 3 the third term represents the volume rate of sediment 
scour (or deposition) per unit length of channel bed. The fourth term 
is envisioned as the time rate of net local exchange between the sus¬ 
pended and bed load zones. An expression defining this term is needed 
for solving Eq. 3. For simplicity, the following parametric exchange 
equation is adapted (Whitham, 1974) 


3 r 

at 


k A : ; 
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(4) 


Here k Q and are constants, T is the concentration at transport capa¬ 
city, corresponding to which, k is the net sediment exchange between the 
suspended and bed zones. Eq. 4 was so constructed in order to (i) 
incorporate in the solution some of the nonlinearity undoubtly present 
in the exchange process, and (ii) preserve the hyperbolic character of 
the sediment continuity equation. Although Eqs . 3 and 4 can be solved 
by successive iterations to obtain c and r, a simpler solution results 
from the following observation. Very near equilibrium the right hand 
side of Eq. 4 vanishes for all practical purposes. That is 


J . 10 




k Q A[(T-c)r-k 1 (R-r)c]^ 0 . 


(5) 


At the same time T deviates slightly from its equilibrium value. Thus, 
from Eq. 5 results 

3r = k l RT 3c . (6) 

at [T+(k 1 -l)c] 2 at 

Also this expression is assumed to hold in slowly varying flow condi¬ 
tions. Approximating Q g by AVC over a small time interval, and combin¬ 
ing Eqs. 3 and 6 gives the sediment continuity equation used in this 
model 



3c „ 3c _ q st 

3t s 8x A(l+f ) 

s 

(7a) 

where 

v = v 

s (1+f s } ’ 

(7b) 


kjRT 

f s " AlT+tkj-Dc] 2 * 

(7c) 


q st = % ' °- A)W If 

(7d) 


Eq. 7a is a quasilinear hyperbolic equation governing the propagation of 
the sediment concentration waves. Its right hand side represents the 
lateral sediment inflow contributed by runoff, tributary channels, and 
bed scouring. The limited experience gained so far with the model 

-3 

indicates that k^ is of order 10 . Thus Eq. 7c has been approximated 

by, 


f £ 
s 


CEL _ , 

AT[1-0.999c/T] 2 


(7e) 


where CEL is a user supplied parameter. This parameter should be ad¬ 
justed to match the time of arrival of the observed sediment load peak 
at the channel outlet. Eqs. 7b and 7e show that the celerity of the 
sediment wave, V , is controlled by the local hydraulic conditions and 
sediment properties since capacity depends on these parameters as well. 
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Eq. 7a can be solved by means of the method of characterist ics. 
From this equation and the total differential of the sediment 
concentration, the following characteristic equations result 


= V 

dt s ’ 
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dt “ A(l+f ) 
s 


equations show that 

Eq. 7a possesses only one 

system of forward 


characteristics. Accordingly, this equation cannot be used in situa¬ 
tions where there are flow reversals. Integrating Eqs. 8a and 8b with 

the initial condition c = c(x ,t ) gives 

o o o 

x 

c = c Q + J (q st /Q) df, (9a) 

X 

o 

x -1 

t = t + / V (1+1 Jdt . (9b) 

o s 

X 

0 

These integrals are used to track the evolution of the sediment waves 
across the characteristic plane. The concentrations existing on all the 
characteristics at the time of their arrival to the downstream boundary 
define the outflow sedimentgraph. 

It is interesting to note that the preceeding equations reflect the 
expected behavior. For instance, it is a recognized fact that the 
streamwise velocity of sediment particles always lags behind the velo¬ 
city of the surrounding fluid (Francis, 1971). This velocity lag ranges 
from very small values for silt particles in suspension, to quite large 
differences in the case of coarse sands and gravels. Thus, the celerity 
of a sediment-1oad wave will be smaller, in general, than the celerity 
ol the carrying flow wave. This trend is also predicted by Eqs. 7b and 
7e, for they indicate that waves of coarse material will travel slower 
than waves of finer material given th.it the carrying capacity of a 
stream increases as the sediment size decreases. 

This celerity lag is strictly a function of local flow and sediment 
properties, and it should not he contused with the differences sometimes 
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observed between the time of arrival of the flow and sediment load 
peaks. In fact, the sediment load may peak ahead, in phase, or behind 
the flow peak depending on antecedent conditions, intensity of the 
event, sediment source location, season of the year, etc. (Guy, 1970). 
To illustrate this point consider a channel reach having a deposit of 
very fine sediment on the upstream portion of the reach, and receiving 
an inflow as shown in Fig. 2a. The sediment will be entrained as 
soon as reaches a sufficient intensity, and will move along charac¬ 
teristics parallel to the flow characteristics (i.e., V g = V, Eq. 8a). 
As the characteristics reach the downstream boundary the water and 
sediment outflows begin to rise, generally at different rates. The 
concentration increase rate is controlled by the ratio of sediment 
entrainment to flow discharge (Eq. 8b). If the supply of loose sediment 
is finite, the sediment outflow will peak at some time t^, while the 
water outflow will continue to increase and peak at a later time t^. 
The lag between these two peaks will obviously depend on the magnitude 
of the sediment supply, inflow rates, and distance of sediment source to 
channel outlet (x = , Fig. 2a). For simplicity, this example has 
ignored backwater effects to restrict the flow movement to forward 
characteristics. Next, consider the same channel reach but with the 
source located farther upstream and a constant base flow 0^ (Fig. 2b). 
The flows and sediment loads reaching the channel inlet (x = x^) are 
generated by characteristics emanating from the upstream reach, x q ^ x ^ 
Xj. Lateral runoff, q^, is assumed to begin entering the channel at t = 
t^. At this point the outlet hydrograph will begin to rise and even¬ 
tually peak at t = t,. . On the other hand, the sediment characteristics, 
which enter the channel reach after t^, will reach the outlet later than 
t^ and will finally peak at a later time t^ > t^. Thus, in this case 
the sediment peak lags significantly behind the flow peak as result of 
changes in the channel water inflow and sediment source location. 

Consider now a channel with sediment moving mostly as bedload. In 
this case the term 9r/3t in Eq. 3 can be neglected. Furthermore, 
Mahmood (197 1 )) has shown that ignoring the time derivative of the spa¬ 
tial concentration has little influence on the simulation of bedload 
propagation. If in addition the period of the flow velocity field is 
very large, the surface gravity waves can be filtered out and only the 
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bedload propagation waves need be retained (Gradowczyk, 1968). In such 
an instance, Eqs. 1, 2, and 3 reduce to the equations governing the 
"known discharge" approximation frequently employed in modeling propaga¬ 
tion of bed transients (Gradowczyk, 1968; de Vries, 1971; Cunge and 
Perdreau, 1973; Thomas and Prashun, 1977; Ponce et al., 1979). 

The foregoing observations demonstrate that the equations presented 
in this Section incorporate the ingredients needed to simulate the 
actual physical process. These equations are complemented with a number 
of ancillary algorithms discussed in the following sections. 

In natural channels there is usually a wide gradation of sediment 
sizes. Although the finer particles comprise the bulk of the load, the 
channel evolution is controlled by the coarser particles (i.e., armor¬ 
ing). Moreover, different sizes are transported at different rates as 
pointed out earlier. It is thus important to predict the movement of 
the individual particle sizes encountered in the channel. The present 
model divides the size range in a suitable number of size fractions, and 
then uses Eqs. 7b, 7d, 7e, and 9 to track the movement of each fraction. 
2.2 ANCILLARY ALGORITHMS 

2.2.1 Sediment Transport Formulas 

These formulas are used to determine the potential carrying 
capacity or a specific flow. Different capacities can be expected for 
different particle sizes, and not all L.isting formulas perform equally 
well for all sizes. In the present model several formulas are used that 
are framed for easy use in digital applications, and require only 
information on the hydraulic parameters of the carrying flow. They are 
the total load formula of Yang (1973) used to estimate the transport of 
very fine to coarse sands (0.1-2 mm), a duBoys-type bedload formula 
(Graf, 1971) used with fine gravel particles (2-4 mm), and the Meyer- 
Peter and Muller bedload formula (Meyer-Peter and Muller, 1948) used 
with particles in the medium to very coarse gravel range (4-65 mm). 
These formulas are presented in Addendum 2 in the forms used in this 
model. The Yang formula was found to give very reliable estimates for 
flows carrying sands in the indicated size range, including flows 
transporting sediment mostly as bed load (Alonso et al., 1980). The 
duBoys-type formula and the Meyer-Peter and Muller formula were selected 
because they were developed from data in the specified size ranges. 


J. 15 





Nevertheless, the user may replace these formulas by others he may deem 
more appropriate for a particular simulation. In particular, simple 
relationships between the sediment transport rate and the flow condition 
developed from in situ field surveys should be given preference. 

2.2.2 Residual Capacity. Composit ion of Material in Transpo rt 

Whenever the bed material consists of a mixture of different 
sediment sizes, the particle size composition of the sediment discharge 
usually differs from that of the bed. Therefore, the transport must be 
characterized by a load rate and a size distribution analysis. Although 
a flow may have the potential to transport sediment, its residual 
capacity or ability to carry any additional Load depends on the sediment 
material already present in the flow. Consider, for instance, a flow 
carrying a load c^ of uniform size d^ and let be the corresponding 
potential capacity of the flow. Then 


T c i = T i - = T , - T i <V T i> 

is the residual capacity of the flow for that size material. The last 
term in this expression represents that portion of 1 already consumed 
by the material in transport. Similarly, if the later were of a dif¬ 
ferent size d 2> the residual capacity for the size d^ materia. , in the 
presence of the load c^, would ne 

T rl = T 1 - T 1 'W ' 

where c 2^2 ma ^ env ^ s i° ne ^ as Chat part of T^ depleted by the load 
c^• If both sizes were simultaneously present in the flow, then 

T rl = T 1 ' VW • T 1 (C 2 /T 2 ) ' 

This expression can be generalized to any size fraction, d^ , and for an 
arbitrary number, n, of load fractions Cj, c^, ...» c , as follows 

T ri = T i - WV - WV 

- T.(c./T.) - •'’ - T.(c n /T n ), (10) 

or, 
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ri 
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i = 1 » 
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( 11 ) 
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n 

The quantity within brackets, 1 - 1 (c./T.), represents the portion of 

j=i J J 

the potential capacity T^ taken up by all the size fractions in 
transport. Therefore, the quantity within brackeLs, A, is the remaining 
capacity for transporting additional material of size d^. For a given 
sediment load, A is the same for all fractions. T. depends uniquely on 
the local flow and the properties of the d^ fraction, while T . depends 
on all these parameters and on the size composition of the sediment 
load. 

When A > 0, any size fraction available for entrainment at the bed 
surface, and for which / 0, will be removed by the flow and added to 
the same sediment size class already in transport. Thus, A > 0 identi¬ 
fies an eroding bed condition. Similarly, when A < 0 the stream carries 
a load in excess of its potential capacity and will deposit the excess 
sediment material on the bed. Therefore, A < 0 characterizes an ag¬ 
grading bed condition. When A = 0 there is no load change and the 
transport process remains in a pseudo-equilibrium condition. By contin¬ 
uously tracking the value of A, the dependence of the individual resi¬ 
dual capacities on the composition of the sediment load and its exchange 
with the bed is readily simulated. 

The size composition of the total load is continuously up-dated by 
adjusting the concentration of individual fractions according to the 
composition of the material removed from or added to the bed. This 
procedure is explained later in the report. 

2.2.3 Bed Compo sition 

In cohesionless beds the material available for entrainment is 


essentially that exposed at the bed surface. As dunes and ripples move 
slowly downstream they continuously mix all the sediment they contain. 
The space occupied by those bed features may thus be regarded as a 
mixing zone below which the bed material remains undisturbed (Fig. 3a). 
Where the sediment contains a mixture of different sizes, the slowly 
moving coarse material tends to collect at the base of the mixing zone 
thus forming a lense of large grains. Under some flow conditions the 
coarse material in the bed may become immobile. In this case, the flow 
washes the finer particles out of the mixing zone leaving an armor co at 
that protects the underlying material. During the scouring of the finer 
materials the exchange between bed and flow takes place in a thin layer 
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Fig. 3. Schematic representation of bed processes 









at the bed surface, identified here as the active layer. Several of 
these layers will be scoured away while the mixing zone is degrading. 
When the bed is armored, the last active layer becomes the armoring 
coat. Conversely, during the process of bed aggradation several active 
layers will be deposited on the bed forming a new mixing zone. 

To model the above processes the mixing zone is pictured as a band 
of constant thickness divided into several layers (Fig. 3b). The layer 
in contact with the flow is always referred to as the active layer. The 
thickness, porosity, and size distribution of this layer can vary 
throughout the simulation, but the layer is assumed to be homogeneous 
within itself at any given time. 

When all the material within the active layer moves, a reasonable 
upper bound of its thickness is obtained, from volumetric considera¬ 
tions, as 


ALT 


100 

P 

n 



( 12 ) 


Here d^, P^, and A^ are the size, percentage, and porosity of the 
coarsest fraction in the sediment mixture. For instance, in a uniform 
bed material with A = 0.50, Eq. 12 gives ALT = 2d^ which is in agree¬ 
ment with the bed load thickness proposed by Einstein (1950). However, 
when the active layer becomes an armor coat, its actual thickness may be 
considerably less than that predicted by Eq. 12 because the bed may be 
armored by particles smaller than d . In that case, it is proposed to 
compute the layer thickness from 


ALT 


100 
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2 P 
i=£ 



1 -A 
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(13) 


where d^ is the smallest grain fraction the flow cannot transport (i.e., 
T.=0, i = £,£+1,...,n). Eq. 13 is also a better measure of the active 
layer thickness when some of the fractions in the bed layer cannot be 
eroded by the flow. At low discharges only the smaller fractions will 
be set in motion and Eq. 13 will thus predict a thin active layer. This 
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is intuitively correct since little bed material will be scoured by a 
low flow during a simulation time step. As the discharge increases, the 
coarser fractions, usually present in small percentages, will be en¬ 
trained and Eq. 13 will predict a thicker layer. This behavior is in 
agreement with the fact that a greater depth of bed can be sorted 
through by a higher flow in the same amount of time. Eq. 13 is thus 
adopted as the general expression for the active layer thickness. This 
equation incorporates Eq. 12, in the limit, by letting S. = n when all 
fractions are in motion. The sediment contained in the active layer is 
the only material available for erosion during a simulation time step. 
When the bed is armored no erosion can occur until the flow develops the 
necessary T for the smallest size fraction present in the armor coat. 
When this happens the armor coat becomes again an eroding active layer. 
If deposition of a certain amount of sediment occurs during simulation, 
this material is added to the bed and a new active layer thickness is 
computed based on the new mixture composition. 

In order to account for the time evolution of the active layer 
thickness, it is necessary to continuously track the size composition of 
this layer. This is done by introducing the following accounting algo¬ 
rithm. Consider a well mixed active layer with three size fractions 
C *l << ^2 << ^3 eroded by a flow with sufficient transport capacity to 
scour all three fractions (Fig. 4). The bed scouring is imagined to 
begin with the entrainment of the d^-particles exposed at the bed sur¬ 
face. Next, the d^-particles at the bed surface are removed, followed 
by the d^-grains hidden underneath the d^-grains. Finally, the d^-par- 
ticles are washed off the bed surface followed by the d^ and d^'parti- 
cles underneath them, and by the d^-grains uncovered by the removal of 
the d^-grains. Thus, one d^-particle is removed for every d^-particle 
entriined by the flow, and one d^ and two d^-grains are associated with 
the removal of every d r ^-grain. This ordering can be readily extended to 
» y number, n, of size fractions, and il is summarized in the following 
"entrainment frequency" matrix. 
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Fig. 4. Schematic reprcspnt.it ion oi active layer composition 
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F = [F..] = 
ij 


2 i-2 2 i-3 2 i-4 2 i-5 2 i-j-l 


where i, j = 1,2,3, ..., n. In this matrix the diagonal elements indi¬ 

cate that every size fraction at the bed surface is scoured once. The 
off-diagonal elements in each row indicate the number of times each 

fraction d., l^i^i, becomes available for entrainment once d. is removed 
J i 

from the bed. Alternatively, the elements of any column, j, express the 

number of times the d^ fraction is depleted due to the entrainment of a 
fraction d^. Obviously, the set of entrainment events in Eq. 14 is one 
of many possible distributions since in actuality any number of d.-par¬ 
ticles can be associated with the removal of a d^-grain. In a more 
general stochastic representation the F matrix would be replaced by the 
conditional probability matrix of the entrainment process. Neverthe¬ 
less, the proposed algorithm is adopted for its deterministic simpli¬ 
city. It should be noted that the F matrix is time invariant and de¬ 
pends only on the number of fractions used in the simulation. 

The amounts of eroded materials corresponding to the frequencies 

F.. depend on the volumes of the fractions contained in the active 
ij 

layer. To convert those frequencies to volumes let be the total 

volume of the d^-size present in the layer per unit channel length, and 

v. . the portion of V. that becomes available for entrainment when the 
ij J 

d^-fraction, occupying the volume , is eroded. Thus, 

v. . = F.. V. , 
ij iJ i 


V. = I F . V , 
J rj r 

j r=J 
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from which, 


v . . 
i-l 


F. . V 
- U 1 


z 

r=j 


F V 
rj r 


x V 

J 


(15) 


From this expression one obtains 
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I v. . 
= i 1 J 
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Eq. 15 can be rewritten as 
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are the percentages of active-layer material in each size class inter 
va 1 . 


Eq. 17 defines the' element:, of a square matrix identified in here 
as the "volume entrainment matrix," v. The order of this matrix is 
equal to the total number of size classes in the active layer. Its 
diagonal elements represent the volumes of the individual fractions on 
the bed surface. The off-diagonal row elements contain the individual 
fractional volumes exposed bv the erosion of larger sizes. Adding up 
these elements gives the vojume of potential erosion associated with the 
removal of the largest size in the row. On the other hand, summing all 
tne elements in each column yields the total volume of each size class 
present in the active layer and available for scour (Eq. 16). These 
concepts are schematically illustrated m Fig. 5a, which depicts the 
v-matrix associated with a small cluster of bed particles grouped in 
three different size rlasses. The shaded area in Fig. 5b represents all 
the material that could be entrained along with tile third fraction. The 
shaded portion of Fig. 5c '(-presents instead the total volume of the 
first fraction contained i n the cluster. During bed degradation, or 
aggradation, the elements . I the- '.-matrix are continuously adjusted as 
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explained in the next two sections. At the end of each simulation time 
step the total volumes of the individual fractions left in the active 
layer are introduced into Eq. 18 to calculate the size distribution of 
the layer. 

2.2.4 Bed Ero sion. Armoring 

Whenever A > 0 the bed is in a degrading mode. This is 
conceptually modeled by depleting the elements of the v-matrix, one row 
at a time, beginning from the smallest fraction. The diagonal elements 
are depleted first, for they are exposed to the flow when simulation 
starts. Next, the off-diagonal elements are scoured one by one 
beginning from the smallest size. When the flow residual capacity is 
not sufficient to remove all the > Junes in a particular row, equal 
amounts (for simplicity) are removed from all the elements in the row 
until the residual capacity is satisfied. The eroded volumes are added 
to the individual fractions already in transport, and the value of A is 
recalculated. This process continues until either A is no longer 
positive definite, or the entire active layer is worked through. For 
example, the volumes (v ), (^v^, ^ V 21 ’ ^ V 22’ ^ V 21 ’ etc -^> ^ v 33^* 
V 31 / 3 ’ v 32 ,/3; V 33^’ v 31 / ^’ v 32 etc -^ would be eroded, in this order, 
out of the v-matrix shown in Fig. 5. However, the extent to which these 
volumes are actually entrained depends on the degree of detachabi1ity of 
the sediment particles. These concepts are used to construct the 
following algorithm expressing the total volume of sediment eroded out 
of each size class in the active layer during a simulation time step: 


2 e.. , e.. = ERO . v.., if AT . ? 2 v. ., (19a) 

j.j iJ ij 


1J ' ri ' j = l 1J 


r< i 

E* =•{ 2 e . . , e . . = <£ 

1 * j=i,l 1J ‘ J 


1 r<i 
tERO.v.., if jSr, AT .=2 v.., 

t iJ ri iJ 


AT . > ‘ v. . , 

r l i U 


(19b) 


r<i 

ERO. AT ., if j = r, AT . = 2 v.., 

r.)’ J r i j = ] ij’ 

AT .< 1 v . , 


( 19c) 



i=l, 2, n. In these equations ERO is a user supplied erodibility 
parameter. This parameter governs the actual amount of bed material 
available for erosion during a simulation time step. ERO is calibrated 
by fitting the sediment yield volume to observed data. Everytime a new 


e_ is computed the concentration of the corresponding load fraction, 
Cj, is upgraded by letting 


. . 

C* = C . + —. (20 ) 

J J A 


This concentration is then entered in Eq. 10 to update A. When A be¬ 
comes nonpositive for any particular size fraction the transport is 
termed "capacity limited." On the other hand, if A remains positive 
after depleting a fraction the situation is termed "supplied limited." 
In particular, if A > 0 after considering all fractions present in the 
active layer, the bed is said to be armored and the active layer becomes 
an armor coat. Fig. 6 summarizes the foregoing simulation sequence. 

After the active layer has been worked through by the flowing 
water, its volumetric composition is given by V^-Ev, i=2, 2+1,...,m, 

where 2 and m are the smallest and largest fractions left in the layer. 
The actual thickness of this material is then 

m V. E* 

AlT * = V * ' (21) 

1=2 i 

Whenever ALT" is equal, or larger, than the thickness obtained from Eq. 

12, the later is taken as the new thickness of the active layer. In 

some instances, however, ALT* turns out smaller than ALT. In these 

cases, a thickness, 6, of undisturbed material is added to the active 

layer (Fig. 7a). Because the model does not track the composition of 

the undisturbed layer, this is always formed by original bed material. 

This material contains, in general, a smaller percentage of the largest 

size class than the ALT* layer (Fig. 7b) and, therefore, 6 must be 

reduced to account for the difference. For simplicity, a simple linear 

correction is assumed given by the ratio between the percentages, and 

P u , in the eroded and undisturbed layers (Fig. 7b). Thus, the corrected 
m 

thickness of the new active layer becomes 
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Fig. 7. Adjustment of active-layer thickness 
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( 22 ) 


P a 

ALT = ALT* + 6 . 

c p u 

m 

The same criteria is used when updating an armor coat thickness, except 

that in this case the ratio P a /P u is replaced by P^/P 1 ?. 

mm r £ £ 

Finally, the local bed elevation of the channel bed is updated by 
subtracting from it the thickness of the eroded material, giving (Fig. 
7a) 


n E. 

z = z - - 1 TV-—-. • (23) 

new old W . , (1-A.) 

i = l i 

2.2.5 S ediment D eposit i on 

When the model senses deposition (A < 0), the flow drops sediment 
on the bed during the time step. At, and the settled material is added 
to the active layer. Within the present deterministic framework, 
deposition begins with the largest sediment fraction and continues 
through the smaller fractions until either the stream is no longer 
overloaded (A = 0), or all the fractions in transport have been 
depleted. The volume deposited out of any fraction, D . cannot exceed 
its (defect) residual capacity, that is, 

Maximum D. = T = T. A . (24) 

i ri i 

However, whether this amount will reach the bed during the time step At 
depends on this being not less than the average time for the sediment 
particles to settle to the channel bed. Data by Jobson and Sayre (1970) 
and by Lean (1971) indicate that this settling time may be computed 
using the particle fall velocities in the quiescent fluid. Therefore, 
the actual deposition of a size fraction during the interval At is 
calculated from 
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particles fall simultaneously from the water surface. The deposited 
volume is subtracted from the material in transport to yield the 

size-class concentration 


From the preceding equations the deposition loop shown in Fig. 8 is 
constructed. 

The settled sediment is added to the active-layer fractions, and 
both materials are assumed to mix thoroughly yielding the volumetric 
composition 


V., when D. = 0, 
1 l 


I V. + D., when D. * 0, 
i i t 

i = 1, 2, . . . , n. These volumes are introduced in Eq. 17 to update the 
bed composition, and this is then used in Eq. 12 to compute the new 
active layer thickness. Finally, the local bed elevation is updated by 
adding to it the thickness of the settled material yielding 
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NUMERICAL SCHEME 


To run the model in a channel reach of length L, the open domain [0 
£ x $. L] x [t £ 0] is covered with a rectangular grid with lines paral¬ 
lel to the x- and t-axes. There is a single family of characteristics 
associated with Eq. 7a. This leads to the solution of an initial, 
one-point boundary value problem. Initial conditions are cross-section 
geometry, bed elevation, average flow velocities, flow depth, active 
width, and bed size composition. Upstream boundary conditions are time 
histories of water and sediment inflow, and size composition of sediment 
load. Let Ax and At be the length and time increments, respectively, 
separating the grid lines (Fig. 9). The coordinates of the grid nodes 

are x = mAx and t = nAt, m, n = 0, 1,2, ... The value of any var- 
m n 

iable, say A, at the node (x , t ) is designated A . Given the above 
’ 3 m’ n ° m 

data specified along the line t = t , the modi is used to compute the 
sediment discharge, load composition, bed elevation, and active layer 
composition at all nodes on the next line t = t n+ |- For use in the 
numerical scheme, the sediment transport equations and ancillary algo¬ 
rithms are implemented as shown in the diagram of Fig. 11. 

The volumes of the sediment fractions available for erosion during 

At are obtained from the active layer thicknesses on the line t = t as 
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v n 

i ,m 


w 11 
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ALT 


100 
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1 » l 

1-A. 


i=l,2, 


n. 


(29) 


These volumes are used in Eq. 17 to compute the volume entrainment 

► • n 
matrix v . 

m 

The flow routing scheme supplied by the user is used to obtain the 
hydraulic-parameter values at the nodes J and K (Fig. 9), and their 
averages are used to compute an average transport capacity, T. The 
concentration of each load fraction at the end of the characteristic 
path is obtained from the following discrete form of Eq. 9a 


n 

c. - c. + 

l ,B l ,m 


q Ax 
^st 


(30) 


where Q is the average of Q T and Q , and q is the piecewise uniform 

J K St 

lateral sediment inflow over Ax. The amount of sediment entrainment, or 
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deposition, and the updated load concent ral i <>ii . c •. are calculated 

using either the erosion loop (Fig. 0j or the deposition loop (Fig. 8) 
depending on the sign of A. If the character. st i i of a particular 
fraction does not pass through the node K (Fig. 9 j, the ioncent rat 1 on at 
this node is obtained by linear interpolation from 


n+1 

C i,m+1 


=< 


... At 
i,B At*' 


Ax 

i, B Ax* 


.., At V 

i , B Ax(1 + f 

s , 1 




Ax(1 + f ) 

.' .s ,J 

At V 


if At* > At, 


i t At* < At. 


( 11a.) 


(31b) 


In these equations Ax* and At* are given by discrete forms of Eq. 9b. A 
similar correction is applied to E* (Eq. 19) to determine the actual 
amount of material eroded in the interval At. The new concentrations 
are used to calculate the new volume load fractions 


n+1 t n+1 

i,m+l C i,m+1 


(32) 


and these are introduced into Eq. 18 to update the load size composi¬ 
tion. The foregoing load updating sequence is summarized in Fig. 10. 
Last, the new bed size composition and elevation are computed. When a 
pseudo-equilibrium condition is encountered (A = 0) no bed updating is 
necessary. In continuous simulation the above calculations are per¬ 
formed consecutively for each channel segment to update all variables 
over the length of the channel. The process is then repeated for each 
time interval in the simulation period. 
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MODEL TESTING 


Tests were conducted to verify the ability of the model to simulate 
various instream processes. Published laboratory and field studies were 
scanned for data suitable for testing the algorithms describing proces¬ 
ses of bed scour, armoring, and unsteady transport. Three useful sets 
of data were found. The first set was collected by Ashida and Michiue 
(1971). They performed a series of laboratory experiments to study the 
effect of sediment gradation on channel armoring. The second data set 
was collected by Lane and Carlson (1933) in the San Luis Valley canals 
in south-central Colorado. They made measurements on the bed and bank 
materials that formed the canals. The beds of these canals have become 
armored over the years. The data offer an excellent opportunity for 
testing the model on a natural system. Finally, the model was tested 
using data collected by the U. S. Geological Survey on a reach of the 
East Fork River, Wyoming (Mahoney et al., 1976). These data offer the 
opportunity of checking the model algorithms on an alluvial stream under 
conditions of unsteady flow. These tests are discussed below. 

As is usual in channels with material consisting of a wide range of 
grain sizes, the bed slopes of the channels used in the above three 
studies were fairly steep. In these cases the kinematic-wave 

approximation to the equations governing unsteady flow of water is 
applicable. In this approximation the momentum equation, Eq. 1, becomes 
Q = KIN . A 1/2 P _1/2 , (39) 

where 


KIN = C 



1-49 r 1/6 
n 



(34) 


is a kinematic-wave parameter. In these relationships P and R are the 
wetted perimeter and hydraulic radius of the channel cross section, 
respectively, C^. is the Chezy coefficient, and n is the Manning 
roughness factor. Eqs. .33 and 34 are also valid when the flow is 
uniform and steady. For the purpose of simulating the East Fork River 
data Eqs. 2, 33, and 34 were solved using the kinematic-wave routing 

scheme developed by Borah et al. (1980). 

The simulations required calibration of the model parameters to 
obtain best-fit of model predictions to observations. One flow routing 
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parameter, KIN, and two sediment transport parameters, CKL and ERO, were 
available for adjustment. CEL controls the travel time of sediment 
load, and ERO governs the total quantity of bed material available for 
entrainment at the bed active layer. 

4.1 FLUME ARMORING STUDY 

Several laboratory experiments were carried out by Ashida and 
Michiue (1971) to study the effects of bed armoring on channel 
degradation. Various sediment mixtures were used with different mean 
diameters and geometric standard deviations. These mixtures were placed 
in a recirculating flume with a bed 2.62 ft. wide and 6S.S ft. long. In 
each experiment a steady uniform flow was passed over the bed. The 
hydraulic conditions were chosen to purposely induce armoring. No 
additional sediment was introduced into the stream, and the rate of 
sediment collection in a trap at the end of the flume was used as a 
measure of the total sediment discharge. 

Data from Ashida and Michiue's run No. 2 were chosen to check the 
performance of the model in simulating the transport of cohesionless 
sediment mixtures and the formation of armoring. The hydraulic condi¬ 
tions used in the run are given in Table 1. In running the model, the 
flume was represented by four reaches of equal length with a point 
inflow of water at the upstream end of the channel. The initial 
particle size distribution of the bed material employed in the run is 
shown in Fig. 12. In simulating sediment transport, the graded bed 
material was represented by ten discrete particle size fractions. They 
are listed in Table 2. Fig. 12 shows the size distributions of the 
eroded material leaving the flume and of the armoring layer as reported 
by Ashida and Michiue, and the best-fit curves predicted by the model. 
The agreement between the simulated and the measured armored layer size 
distributions is good. However, some discrepancy exists between the 
curves for the eroded material. The measured curve definitely shows 
larger sizes present in the sediment load. This discrepancy may be 
attributed to inaccuracies in the estimation of sediment entrainment. 
In the transport capacity formulas an average tractive force concept is 
implied, which excludes the effect of the large instantaneous 
fluctuations associated with turbulent tractive forces (Alonso and 
Coleman, 1981). 
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Table 1 - Hydraulic conditions in the flume and San Luis Valley Canal 
studies. 


Data 

Source 

Flow 

Rate 

__ (cfs) 

Flow 
Depth 
(ft ) 

Average 

Velocit y 
(f£s) 

F.nergy 

SI ope 

Manning's 
Coefficient 

Ashida and 
Michiue (1971) 
Run 2 

1.06 

0.22 

1.81 

0.00440 

0.017 

Lane and 

Carlson (1953) 

128.0 

1 . 77 

4.00 

0.00243 

0.023 


Test Section 12 


Table 2 - Size distributions used for simulating the flume and San Luis 
Valley Canal data. 



Flume 

armoring study 

San Luis Valley 

Canal 

Size Class 

Percent 

Size class 

Percent 

Interval 

per class 

Interval 

per class 


(mm) 

interva1 

(mm) 

interval 

0.2 

- 0.3 

9.5 

0.000 - 0.149 

4.0 

0.3 

- 0.4 

9.5 

0. 149 - 0.297 

5.7 

0.4 

- 0.6 

11.0 

0.297 - 0.590 

10.1 

0.6 

- 0.8 

6.0 

0 590 - 1.190 

12.0 

0.8 

- 1.0 

4.0 

1.190 - >.380 

9.2 

1.0 

- 2.0 

12.0 

2.380 - 4.760 

5.4 

2.0 

- 4.0 

18.0 

4.760 - 9.525 

13.3 

4.0 

- 6.0 

18.0 

9.525 - 19.05 

17.9 

6.0 

- 8.0 

10.0 

19.05 - 38.10 

13.2 

8.0 

- 10.0 

2.0 

38.10 - 76.20 

9.2 
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fig. 12. Size distribution ci.r/is oi'i.lined for f J nine armoring test 









4.2 


SAN LUIS VALLEY CANAL TESTS 


The test described in this section is based on data collected by 
Lane and Carlson (1953) on a reach of one of the canals located in the 
San Luis Valley of southern Colorado. These canals were constructed in 
the late 1800's in the alluvial cone deposited by the Rio Grande River 
on the floor of the San Luis Valley. Near the appex of the cone the 
deposits consist of sand, gravel, and cobbles, the size of the cobbles 
decreasing with the distance from the appex. The canals were found very 
stable, and therefore presented an unusually favorable condition for 
studies of stable canals. Several reaches ranging in length from 600 to 
2200 ft. were selected for study. Hydraulic measurements were made on 
the most regular portions of the reaches. Observations and mechanical 
analysis were also made of the bed and bank materials forming the test 
sections. Samples of the material in which the canal was constructed 
were obtained by excavating into the bank of each section. Mechanical 
analysis of the material forming the bed surface layer disclosed that in 
the stable sections the finer material had been removed from the layer 
and an armor coat of coarser material was left. 

Data from the test section No. 12 was selected to simulate the 
development of the armor coat. The hydraulic conditions observed in 
this section are presented in Table 1. The bank and bed material size 
distributions are shown in Fig. 13. For simulation purposes, a 600 ft. 
long reach divided into four equal segments was assumed, with a constant 
inflow of clear water. The initial bed-material size composition, 
assumed equal to that of the banks, was represented by the distribution 
listed in Table 2. The size composition of the active layer at 2.5, 
12.5, 37.5, and 87.5 hours are shown in Fig. 13. As time increases, the 
layer is quite rapidly depleted of particles 10 mm in size or smaller, 
and the distribution in this size range approaches assymptotically the 
measured curve. However, the simulated distribution in the 10 to 75 mm 
r. i ze range differs from the composition. This difference may be 
explained in part by the use of a constant rate of discharge. Although 
information related to the actual history of flows in the canal is not 
available, Lane and Carlson (1953) mentioned that the canal had 
sustained flows considerably above those observed dining their study. 
Therefore, the observed bed-layer composition was most probably shaped 
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Fig. 13. Evolution of the mzp distribution in to[> layer of San Luis Valley 
Cana 1. 
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by flows 

having transport 

capacit i es 

larger than 
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this 

simulation 

Another reason 
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use ol average 
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section. 

Nevertheless, 
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clear 
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the measured 
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distribution over the entire size range. 

4.3 EAST FORK RIVER PROJECT 

A detailed program of hydraulic and sediment transport data 
collection was undertaken by the U.S. Geological Survey during 1975 on a 
reach of the East Fork River near Boulder, Wyoming. The study reach is 
approximately 3 miles in length and has a tributary area of about 180 
square miles. About half of this area lies within the Wind River 
Mountains. The other half is provided by the area drained by the Muddy 
Creek tributary. A map of the study reach showing the position of the 
principal gaging stations is given in Fig. 14. The river in the study 
reach is about 60 ft. wide and 4 ft. deep at bank full stage. Most of 
the water during high flows comes from melting snow of the mountain 
area. In the East Fork River the bed is gravel on the riffles and bars, 
but coarse sand constitutes the bulk of the bed load. The bed material 
in the Muddy Creek is sand but much finer than the sand in the East 
Fork. Stage recording gages were installed at sections B-l, B-5M, and 
B— 17. Sediment inflow to the reach was measured at sections B-l and 
B-5M by sampling suspended sediment with a standard DH-48 hand sampler 
and determining the bedload discharge with a Helley-Smith bedload 
sampler. The bedload discharge past the section B-17 was measured with 
both a bedload trap and a Helley-Smith sampler. Leopold and Emmett 
(1976) and Mahoney et al. (1976) describe in detail the stream, the 
bedload trap, and the data <ollected during 1975. The data included 
cross-sectional surveys, flow rating curves, water temperature, sediment 
transport rates, and particle-size distributions of bed load material. 
This data set was also simulated by Bennett and Nordin (1977). Their 
input data reduction procedures were followed to some extent in the 
present test. 

In this simulation, the study reach was divided into the fourteen 
subreaches shown in Fig. 14. Surveys of cross sections were used to 
obtain relationships between area, wetted perimeter, and stage. Time- 
averaged active-bed widths were determined from the surveys by comparing 
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plots of bed eLevation at the beginning and end ot the simulation per¬ 
iod, and observing which parts of the bed moved vertically. As noted hv 
Bennett and Nordin (1977) this prodecure may overestimate the active 
width because not always all parts of the active width move simultan¬ 
eously. However, there may be no better way, short of continuous bed 
monitoring. The profile of six of the cross sections, as measured on 
June 2, 1975, are plotted ia Fig. 15. Also plotted are the estimated 
active-bed widths, and the stages observed during June 7. As these 
figures reveal, overbank flow occurred at several sections along the 
reach during the high flow periods, including the inlet sections H-1 and 
B-5M. 

Flov input to the system was provided by point loads at sections 
B-l and B-5M. These point loads were generated by converting ' ! .e stage 
readings to flow rates using the rating curves constructed from the 
stages and flow discharges measured at those sections The rating 
curves reported by Mahoney et al. (1976) are shown in Figs, lba and 16b. 
The sum of the inflow hydrographs measured during the simulation period 
is about equal to the outflow hydrograph observed at section B-17, 
indicating that water was carried by essentially translatory waves. 
This fact lends further support to the use of a kinematic-wave routing 
scheme. 

Inflow of sediment at sections B-l and B-5M were obtained from 
sediment-flow rating curves. These were constructed by relating the 
measured sediment loads to water discharges obtained from the above flow 
rating curves. Suspended sediment did not contribute signii"cantly to 
the total input load and, therefore, the sediment discharges were based 
solely on the bedload data. When plotted in log-log scales the data 
exhibited considerable scatter but displayed a linear trend. Linear 
interpolation of the logarithmic values yielded the rating <urves shown 
in Figs. 16c and l6d. The data scatter resulted in the low correlation 
< ■•fficients indicated in the figures. A total ot ten different frac¬ 
tions with mean diameters 0.12, 0.2 5. 0.5. 1.0, 2.0, 4.0, 8.0, 16.0, 
32.0, and 64.0 mm were selected t <> repr< rut the sediment in the bed and 
in transport. The percentages >1 inuleii.il associated with each size 
were determined from the sampled i/r d ■ ;,t r i but i ■••ns These percentages 
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were used to convert the total bedloads to sediment transport rates for 
each particle-size fraction. Similar analysis was performed for the 
data collected at section B-17. These data are used in comparing 
measured and predicted water and sediment discharge 1rom the study area. 

A thirty minute time step was used in modeling the East Fork River 
for a period of twenty-two days 1roin May 29 through June 19, 197S. The 
first step involved calibration of the flow parameter to best fit the 
observed outflow hydrograph. A comparison of simulated and recorded 
flows at section B-17 is presented in Fig. 17. Agreement is satisfac¬ 
tory at low and intermediate fI<ws, but the measured values exceed the 
simulated values at high flows. The reason for this discrepancy is the 
occurrence of overbank flows which were not accounted for in the con¬ 
struction of the flow rating curves (Figs. 16a and 16 b) . 

The next step after obtaining the hydraulic results was simulation 
of the sediment transport. The simulation was performed a number of 
times, adjusting each time the sediment parameters to obtain the best 
possible agreement between simulated and observed values. Fig. 18 shows 
predicted and recorded bedload discharges at section B-17 for the 
twenty-two day simulation period. The dashed lines joining the data 
points are imaginary lines used to approximate the shape of the measured 
sedimentgraph. The observed total bedload outflow between June 1 and 
June 19 is 1889 tons. The predicted value is 1942 tons, less than three 
percent higher. In spite of this agreement, there is a tendency for the 
simulated bedload rates to be lower than the recorded values at high 
flows, and higher during the recessions. The differences are most 
noticeable during the first high flow period. A possible explanation 
for these discrepancies is provided by field observation in the East 
Fork River recently reported by Meade et al. (1981). Over the years 
they have observed that, during the runoff season, some sections of the 
reach are scoured while other sections are filled, resulting in a 
normniform storage of movable bed material along the stream. 
Consequently, the relation of bedload rate to water discharge varies 
significantly from one part of the reach to another. Meade et al. 
(1981) noted that iminediate 1v downstream of a storage area the bedload 
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transport may increase steeply with the initial increase in water 
discharge and later, as the stored movable material becomes depleted the 
transport rate rapidly decreases, resulting in a multi-valued rating 
curve. They also mention that this behavior has been indeed observed at 
the gaging station B — 17 (Leopold and Emmett, 1976). This behavior 
cannot be exactly reproduced by the present model, which assumes a 
supply of transportable material uniformly distributed along the 
streambed, and single-valued relations between bedload transport and 
water discharge at all sections. Nevertheless, the predicted 
sedimentgraph displays the rising and falling trends observed in the 
recorded bedload during the twenty-two day Lime span. 

Fig. 19 shows a comparison of simulated against sampled size dis¬ 
tribution of the bedload at the beginning, middle, and end of the first 
high flow period. The model predicts fairly well the distribution of 
the coarser materials, and the shift of the d during the event. The 
recorded values, however, indicate less mobility of material in the 
finer size range than predicted. There could be several reasons for 
this disagreement. For one thing the kinematic-wave approximation used 
in routing the flow may be distorting the local energy slopes, resulting 
in overestimation of entrainment threshold conditions. Or, there might 
have been more fine material being carried in suspension which could not 
be sampled as part of the bedload, but which was actually accounted for 
by Yang's total load formula used in the simulation. A definite answer 
requires further investigation. 
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5 CONCLUSIONS AND RECOMMENDATIONS 

5.1 CONCLUSIONS 

1. A one-dimensionaI numerical model has been developed for simulating 
the movement of well graded sediments through a stream network. 

2. Hydraulic routing is performed by using any acceptable algorithm 
supplied by the user because the water movement is assumed 
uncoupled from the sediment process. The model can be used in 
conjunction with any suitable sediment yield model to supply the 
water and sediment runoff from lateral areas. 

3. The sediment routing scheme is based on the physical processes 
governing the mechanics of sediment movement in alluvial channels. 
The model recognizes the effect of bed and suspended load 
interaction on the total load movement, can simulate bed armoring, 
changes in bed elevation, and longitudinal sorting of eroded 
material. 

4. The applicability of the model is restricted to noncohesive 
materials, relatively stable channel geometries, streams with 
negligible in and out-of-bank transport, and flows in which 
transverse currents may be ignored. 

5. The model gave satisfactory results when tested on laboratory data 
from a flume armoring study, and field data from the San Luis 
Valley Canal, Colorado, and the East Fork River, Wyoming. These 
test . tend to indicate that the model adequately simulates the 
transport of graded cohesionless sediments, including the effect of 
armoring. 

5.2 RECOMMENDATIONS 

i . It is recommended that the channel model he further tested against 

a variety of real situations with special emphasis on the scour, 
deposition, and transport of noncohesive materials. 

2. it is recommended that the model he further developed and refined 

to include the following capabilities: (a) improve the one¬ 
dimensional representation by separating flows in the incised 

channel from flows over flood plains; (b) account for in and out- 
of-bank transport, and lateral distribution of bed-material 
properties and hydraulic conditions; (c) predict the variation of 
lateral bed slope and lateral sediment sorting around channel 
bends; and (d) simulate sediment retention by grasses and other 
vegetative covers. 


J . 52 







3. The channel model should be tested using hypothetical situations to 
confirm that the model does respond in a realistic manner. for 
instance, these tests may include the following channel-stabi 1 ity 
related applications: (a) Consolidate the channel model with 
continuous sediment yield and bank-stability models. Run the 
consolidated models for a period of a few years, and predict the 
size and grade of channel needed to maintain a bank height-slope 
that is stable for a given stratigraphic condition. (b) Run the 
consolidated model for a combination of unstable bank and steep 
grade and observe what combination of bed armoring and/or grade- 
control structures are predicted to stabilize the channel. (c) 
Select a range of storm events and use the consolidated model to 
study slough of bank material, and find channel width and/or 
armoring coat that is needed to prevent erosion of slough material. 

4. It is recommended that data gathering efforts be continued to 
provide an adequate base lor further model development and 
va1idation. 
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ADDENDUM 1: SEDIMENT CONTINUITY EQUATION 


Consider an unsteady streamflow carrying sediment down a channel 
with arbitrary cross-sectional geometry and alignment (Fig. 1.1). The 
volume concentration of the sediment-laden flow, c, may vary from point 
to point, and as a result of convection, from instant to instant at any 
point. The lateral volume rate of sediment inflow, q^, can vary with 
both space and time. The continuity equation for an infinitesimal unit 
volume in the neighborhood of a fixed point is 


3c 

3t 


+ 



= o. 


( 1 . 1 ) 


in which u^ is the instantaneous local velocity component of the 
sediment-laden fluid along the x^-direction. Repetition of the sub¬ 
script k in a term implies summation over the three orthogonal coordi¬ 
nate directions. In a turbulent flow, time averaging Eq. 1.1 yields 


3 c 3 

3t 3x. 

k 


(c u k + c’u^ ) = 0, 


( 1 - 2 ) 


where the overbars indicate temporal means, and the primed terms repre¬ 
sent turbulent fluctuations. Let assume the mean cross sectional area, 
A, of the channel divided in two parts. One part, A^, is occupied by 
the sediment carried in suspension, and another, A^, is occupied by the 
sediment transported as bedload. It should be noted herein that the 
orientation of the coordinate system is arbitrary and, in general, A is 
a function of as well as of time. For convenience, the coordinate 
system will be chosen so that the direction normal to A coincides with 
the streamwise direction x (Fig. 1.1). The continuity equation for the 
suspended-load section is obtained by integrating Eq. 1.2 over A^, to 
obtain 


;; 


3c 

at 


- + ff (CjU + rju' ) dA 


‘l 
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Using Leibnitz’s rule: 


r) — — — 1 ^ _ _ _ 

JJ d Aj - [c + 3 X //(CjU + cju’J d Aj 


A1 
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A1 


3A. 


I(c l u + C 1 U ’ J 9x~t = °> 


(1.4) 


in which is the mean perimeter bounding the mean area A^. 
dary condition for this cross section is 


The boun- 


dA 


8A 


9A, 


[c i dT ] = lc i ( aT + 11 3x I)] = - 1 ' c i % d V 


(l.s) 


°i °i 


where q^, the volume rate of lateral sediment inflow per unit length of 
, is taken positive for inflow. Expanding Eq. 1.5 and time averaging 
gives 


3A j _ 9A __ 3A 

[ "i aF + "i “ ax + w ' 9x 


X ( c ! + cjqp do ] . (1-6) 


°) °1 


Substituting this equation into Eq. 1.4 yields 


\l SI CjdAj* SJ (CjU + cju , )dA ] = / (c j q^ + rjqpdOj (1.7) 


1 


1 


1 


Analogous to the Reynolds closure scheme for turbulent momentum trans¬ 
fer, an eddy mass diffusivity tensor, e, is introduced such that 


c:u = -e 


3x 


( 1 . 8 ) 


Now let 


C] - + c», 


(1.9) 
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and 


u = u + u*, 


( 1 . 10 ) 


where the tildes denote spatial averages over A^ , and the asterisks 
denote deviation from the spatial averages. Substitution of Eqs. 1.8, 
1.9, and 1.10 into Eq. 1.7 gives 


3t 3x 


dx 


ssu 



dx 


c*u*)dA + 


where 

Qj = dA 2 

A1 


[ (c l q £ + c l q £ ) d V 

ol 


( 1 . 11 ) 


represents the volume rate of suspended sediment discharge. The first 
integral on the right hand side of Eq. 1.11 denotes the rate of longitu¬ 
dinal dispersion of sediment. The second integral can be replaced by 

the sum of the sediment inflow from runoff and tributaries, q , and the 

*s 

rate of sediment volume transfer across H-H' (Fig. 1.1) from the bedload 
zone, denoted herein as - Or^/Ot. Longitudinal dispersion in the sus¬ 
pended zone is negligible with respect to vertical dispersion and can be 
omitted. Thus, Eq. 1.11 becomes 


9 V~1 

8Qj 

3r 

3t 

+ ~— - q 

ox s 

3t 


( 1 . 12 ) 


Similarly, integrating Eq. 1.2 over A^ yields the continuity equation 
for the bedload, 


3 Vi *2 

3t + 3x 


J (? 2 q 2 + e'q') do 2 , 


(1 -13) 
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where is the volume rate of bedload discharge, and is the mean 
perimeter bounding . Two different sources contribute to the integral 
in the right hand side of Eq. 1.13. One is the rate of sediment volume 
transfer across H-H' from the suspended zone, and denoted Sr^/St. The 
other contribution comes from the rate of sediment exchange with the bed 
surface, which can be approximated by 

q b = -(l-A)W || , (1.14) 


where A is the bed porosity, W defines the active bed width (Fig. 1.1), 
and z is the local bed elevation above a reference datum. The negative 
sign in Eq. 1.14 accounts for the fact that the sediment flux across 
is negative when sediment settles out of the bedload zone during bed 
aggradation (3z/3t>0). Replacing the integral in Eq. 1.13 by the sum of 
S^/Bt and q^ gives 


3A c 3Q 3r . 

-.- 2 2 + __ 2 = — 2 - ( i - X)W — 

3t 3x 3t U Al 3t 


(1-lb) 


Adding Eqs. 1.12 and 1.15 results in 


3Ac 9Q s n 3z f 3r l 

3t~ + 3x~ + (1 ' X)W 3t + [ 3T 


3r„ 

_ i 

3t 


= 


(1.16) 


where Ac = A^Cj + aua Q s = Qj + ^2' Th e tern) within brackets 
represents the net rate of sediment volume transfer across H-H'. This 
term is different from zero whenever a net amount of sediment passes 
from the bedload to the suspended load zone, or vice versa, as a result 
of bed scour or sediment deposition. Letting 


3r 

3t 




and dropping the tildes and overbars, Eq. 1.16 becomes 


3Ac ( , , 3z 3r _ 

9t~ + 3x + (1 * A) W 3l + 3t ~ V 


(1.17) 
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in which Q g is the total volume sediment discharge, and c is the total 
average volume concentration. Eq. 1.17 is the sediment continuity 
equation used in the present model, and equally applies to streamflows 
carrying sediment mostly in suspension or as bedload. 
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ADDENDUM 2: SUMMARY OF TRANSPORT FORMULAS 
This Addendum describes briefly the transport formulas used in the 
present model. Each formula is presented as the originator intended, 
but rewritten in terms of dimensionless parameters. Where the formula¬ 
tions required graphical solutions (i.e., determination of threshold 
conditions from Shield's curve), analytical equivalents (not shown here) 
have been worked out to facilitate their use in digital computation. 
Total load formula of Yang (1973): 

* = ej* 0 Z 50 (V/n*)[l0* * 6 /S), (2.1) 

where 


$ = 5.435 - 0.286 log (w d^Q/v)-0.457 log (u^./ w ) + 
[1.799 - 0.409 log (w d 5 /v) - 0.314 log (u*/w)] 
log (VS Q /w - V c S Q /w), 


v J 2.5/[log(u*d 5Q /v) - 0.06] + 0.66, 0<u J .d 5Q /v<70, 

^ 2.05, u*d 50 /v > 70. 

DuBoys-type bedload formula (Graft, 1971): 


( 2 . 2 ) 

(2.3) 

(2.4) 


where 


<t> = a e 2 (i - e/e ) 


a- x (S-l) 3 / Vg 3 / 2 d 1/2 


(2.5) 

( 2 . 6 ) 


Empirical relationships for the sediment coefficient x ate given by 
Graft (1971). 

Bedload formula of Meyer-Peter and Muller (1948): 

<}> = 8 e 3 £ 2 [k 3/2 - e c /e m ] 3/2 , (2.7) 

in which 

K = (V/u*)(f b /8) 1/2 (2.8) 

The friction factor associated with bed skin friction, f^, is obtained 
from (Vanoni, 1975) 
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37.40R 
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, u k 

H 1 * 3< ^= 7 °* 


u.,_k 


(f b ) ’ = 2 log (4R b /k s ) +1.14, —* > 70, 


(2.9) 

( 2 . 10 ) 


where k g = dg^, Np is the flow Reynolds number, and R b is the bed hy¬ 
draulic radius. 

In Eqs. 2.1, 2.5, and 2.7 4> is the dimensionless volume transport 
rate, 0^ and are the mobility number and relative roughness based on 
the d^ grain size, and the subscript c denotes threshold conditions. 
These parameters are defined as 



* = CQ S /W)/[(S-1)gd| Q ]^, 

(2. 

11) 


e k = 4/[S-l)gd k ], 

(2. 

12) 

and 





>> 

II 

CS3 

(2. 

13) 

where S 

is the specific gravity of the bed material, and u, 

j s the 

bed 


shear velocity. 
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ADDENDUM 3: DESCRIPTION OF COMPUTER PROGRAM 

The computer program used in the simulation of the East Fork River 
data is listed below. The program consists of a main program and 
several subroutines. The main program inputs the required data to the 
system, calls the subroutines according to the computational scheme, and 
prints out the calculated results. Subroutine WROUT is a user supplied 
program that routes water through a channel reach for each time step. 
Subroutine SROUT performs sediment routing through a reach for each time 
step, and it calls in turn subroutines DUBOYS, MEYER, SETVEL, SHIELD, 
and YANG. These subroutines compute potential carrying capacities and 
sediment transport parameters. 

In order to minimize the need for core memory, the program was 
organized by taking advantage of the absence of backwater conditions in 
the present simulations. In the adopted organization the length of the 
entire channel reach is divided into a few segments, or channel blocks, 
and the time steps used for the whole simulation period are grouped into 
several consecutive time blocks. Then, the channel blocks are processed 
in sequence for each time block. Input data for each channel and time 
block are stored in disk files. A list of the important variables in 
the computer program is given in the following section. 

The codes requires 72,741 words on a Mod Comp Classic computer 
system. This is a 16-bit machine that uses two words for each single 
precision variable. The execution time for the East Fork River test is 
approximately 9 minutes. 

List of Fortran Variables 

Name Description Units 

ACCM(IF, IFA) Element of volume entrainment matrix ft 3 

at the start of the current time step. 

IFA indicates the material fraction 
exposed by the removal of fraction IF. 

AE Area of flow cross section at the ft 2 

start of the current time step. 
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ARMHT 

BEDELV(IC) 

BEDMAT(IF) 


BEDUP 


BEDWN 


CAP 


CC(IC, IF) 


CDEP(IC) 


CHEZY 

CHI 

CICIC, IF) 


Current thickness of active layer. 

Bed elevation of section IC at 
the end of the current time step. 
Vector used to store the volume of 
individual material fractions present 
in the bed layer. 

Change in bed elevation caused by 
deposition during the current time 
step. 

Change in bed elevation caused by 
erosion during the current time step. 
Volume concentration of individual 
material fractions at transport 
capacity. 

Volume concentration of material 
fraction IF passing through the 
section IC at the end of the current 
time step. 

Coefficient of power formula relating 
area of cross-section IC to water 
depth. 

Chezy's roughness coefficient. 
Coefficient for the DuBoys sediment 
transport formula. 

Volume concentration of material 
fraction IF passing through the 
section IC at the start of the 
current time step. 


ft 

ft 


ft 3 


ft 


ft 


ft 3 /ft 3 


ft 3 /ft 3 


L 

ft V sec 


ft 5 - 25 / 


ft 3 /ft 3 
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COEF(J) 


Coefficient of the entrainment 


CONC 

CPER(IC) 

CUP(1T, 

DARM 

DELT 

DELX 

DEPO 

DMM(IF) 

DPTH 

DTS 

EDEP(IC) 


IF) 


frequency matrix. 

Capacity concentration predicteii by 
transport formulas. 

Coefficient of power formula relating 
wetted perimeter of section IC to 
flow cross-sectional area. 

Volume concentration of inflow of 
material fraction IF at the start 
of the current time step. 

Size of smallest bed material fraction 
that becomes immobile in the active 
layer. 

Time taken by a sediment characteris¬ 
tic to travel the distance DELX. 
Channel length increment used in the 
computational grid. 

Volume of fraction IF deposited on 
the bed during the current time 
step. 

Representative size of sediment 
fraction IF. 

Water depth. 

Size of current time step. 

Exponent ol power formula relating 
area of cross-section IC to water 
depth 


ppm 


ft 3 /ft 3 


mm 


sec 


ft 


ft 3 


mm 


ft 

sec 
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EPER(IC) 


Exponent of power formula relating 


ERS(IF) 


ERO 


G(IT, I, IF) 


GAMA 

GMES(IT, I, IF) 


GTOTC(IT) 


GTOTM(IT) 


I 

IF 


INLS 


wetted perimeter of cross-section IC 
to flow cross-sectional area. 

Volume of sediment fraction IF 
eroded from the bed during the 
current time step. 

Parameter controlling detachment of 
bed material. 

Calculated volume discharge of sedi¬ 
ment fraction IF passing through the 
downstream end of channel block I, 
at the end of the current time step IT. 
Specific weight of water. 

Measured inflow of sediment fraction 
IF to channel block I, at the start 
of the current time step IT. 

Calculated sediment discharge passing 
through the channel outlet at the 
end of the current time step IT. 
Measured sediment discharge passing 
through the channel outlet at the 
start of the current time step IT. 

Index identifying channel block. 

Index identifying sediment-size 
fraction. 

Number of channel subreaches in a 
channel block. 


ft 3 


ft 3 /sec 


lbs/ft 3 

lbs/sec 


lbs/sec 


lbs/sec 
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IT 


Computation time step. 


ITCOM 

KIN 

GLAT(IT) 


MEYERP 


N 


NARM 


NFR 

NSEG 

PCBC(IC, IF) 


PCF(IC, IF) 


PCFF(I, IC, IF) 


Number of time steps in simulation 
period. 

Kinematic-wave routing parameter. 

Lateral volume inflow of sediment ft 3 /sec 

to channel block, during the current 
time step IT. 

Coefficient in the Meyer-Peter and 
Muller sediment transport formula. 

Number of time blocks the simula¬ 
tion period is divided into. 

Integer identifying the smallest 
size fraction that becomes immobile 
in the active layer. 

Number of representative size 
fractions used in the simulation. 

Number of channel blocks. 

Percentage of material fraction IF 
present in the active bed layer of 
cross section IC. 

Percent of material finer than size 
IF present in the active bed layer 
of cross section IC. 

Initial percent of material finer 
than size IF present in the active 
bed layer of cross section IC in 
the channel block 1. 
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PCW(IF) 


Percent of material in transport 
finer than size IF passing through 
the channel outlet. 


POR(IF) 

Porosity of bed material fraction IF. 

ft 3 /ft 3 

Q(IT) 

Computed water discharge at the end 

of the channel block I, and at the 

end of the current time step IT. 

ft 3 /sec 

QLAT(IT) 

Lateral water inflow to channel block, 

during the current time step IT. 

ft 3 /sec 

QMES(IT, I) 

Measured water inflow to channel 

block I, at the start of the 

current time step IT. 

ft 3 /sec 

QUP(IT) 

Upstream water inflow to channel 

block, at the start of the 

current time step IT. 

ft 3 /sec 

RESCAP 

Residual transport capacity of an 

individual material fraction, 

ft 3 /ft 

- 

expressed as volume of dry sediment 

per unit length of channel. 


RHB 

Hydraulic radius. 

ft 

SCAP(IF) 

Potential transport capacity of 

material fraction IF, expressed as 

volume of dry sediment per unit 

length of channel. 

ft 3 /ft 

SLN 

Length of channel block. 

ft 

SLOPE(IC) 

Channel bed slope at section IC. 

ft/ft 
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SHU 

Average value of the kinematic 

viscosity of water for the simu¬ 
lation period. 

ft 2 /sec 

SPGR 

Specific gravity of sediment material. 

- 

SUM 

Summation term in Eq. 11. 

- 

TAO 

Average unit tractive force. 

lbs/ft 2 

TBM 

Total volume of bed materia] 

contained in the active layer per 

unit length of channel. 

ft 3 

TC 

Critical unit tractive force. 

lbs/ft 2 

TCA(IF) 

Sum of all the elements in row IF 

of the volume entrainment matrix. 

ft 3 

TEMP 

Average water temperature for the 

simulation period. 

Fahrenhei 

TGIN(IF) 

Vector used to store the measured 

volumes of all fractions that enter 

the channel during the simulation 

period. 

ft 3 

TGM(IF) 

Vector used to store the measured 

volumes of all fractions that leave 

the channel during the simulation 

period. 

ft 3 

TGMES 

Total sediment yield measured at the 

channe1 outlet. 

lbs 

TGO(IF) 

Vector used to store the computed 

yield of individual fractions. 

lbs 


J. 71 




TGOUT 


Total sediment yield computed at 
the channel outlet. 


lbs 


UST 

VEL 

VS(IF) 

WATMAT(IF) 

WEP 

WIDTH(IC) 

XSI(IC) 


Bed shear velocity. 

Average velocity of flow. 

Settling velocity in quiescent water 
for material fraction IF. 

Vector used to store the volumes of 
all the material fractions being 
carried by the flow during the 
current time step. 

Wetted perimeter of channel cross 
section. 

Vector used to store the active bed 
width of all channel cross sections. 
Distance of cross-section IC to 
upstream boundary of channel block 
containing IC. 


ft/sec 
ft/sec 
ft/sec 

ft 3 


ft 

ft 

ft 


1 






o n 


C LIST OF COMPUTER PROGRAM USED IN THE EAST FORK RIVER TEST 
C fr******************************************************** 

C 

DIMENSION TITLE(20),QMES(264,3),QOUT(30 0 >,GMES(264,10,3), 

4PCFF <3,6,10),PCF <21,10),PCFR(21,10),GTOTM < 264),CTQTC(264), 

&QS(21),INT(21),GS(21),CII<3,6,10),BEDLVL(3,6),TGO(10 >,TGM(10) , 
6TGIN(10) 

COMMON /ROUT/ I,IT,SLN,CHEZY,DTS,ITCOM,INL,QBASE, 

&QUP(300),QI(S0,2),XI(SO),KSI(SO) ,0(300,3) 

COMMON /WROUT/ SLOP,QLAT(300),QC(SO,2),XC(SO>,KSC<SO) 

COMMON /SROUT/ SPGR,CAMA,SNU,NFR,INLS,WIDTH(10),SLOPE(10 ) , 

6CPER(10),EPER(iO),CDEP(10),EDEP(10),GLAT(300,10),XSI(10), 

&XSC<10,10),BEDELV(10),CI(10,10>,CC(10,10),PCBI(10,10),PCBC(10,10 > , 
6CUP(300,10),DMM(10) ,COEF(10),US(10),POR <10),G(264,3,10 >,PCW<10), 
6ER0,CHI,CEL 
C 

C DATA INPUT 
C GENERAL INFORMATION 

READ(4,406) NSEG,DTS,ITCOM,NFR 
C PHYSICAL PROPERTIES 

READ(4,407) TEMP,GAMA,SPGR ,SNU 
READ(4,408) (DMM(IF>,IF = 1,NFR ) 

C WATER AND SEDIMENT ROUTING PARAMETERS 
READ <4,40V)ERO,CHI,CEL,CHEZY 
C INITIAL BED MATERIAL SIZE DISTRIBUTION 

READ(4,40S) ((PCFF(1,IC,IF> , IF - 1,1 0 > , I C= 1,3 > 

READ(4,405) ((PCFF<2,IC,IF),IF=1,10),IC=1,4) 

READ(4,405) ((PCFF(3,IC,IF),IF-1,10),IC-1,6) 

C COEFFICIENTS OF ENTRAINMENT FREQUENCY MATRIX 
COEF(1)=1.0 
COEF(2)=1.0 
DO 280 J-3,NFR 
280 C0EF(J>-2.OfcCOEF < J-1 ) 

SUBW= (SPGR-1 . 0 ) i'GAMA 
C POROSITY AND SETTLING VELOCITY 
DO 88 1=1,NFR 

POR <I) = 1.-0.245-0.0864/(0.S*DMM(I>>**0.21 
D=DMM(I) 

CALL SETVEL(D,W) 

VS< I)=U/30.48 
88 CONTINUE 

BEGINNING OF TIME-BLOCK LOOP 
DO 999 N=i,4 

C MEASURED WATER DISCHARGE AT SECTIONS B-l, B-5M, AND B-17 
DO 351 1=1,3 

READ(3,333) < QMES(IT,1),IT = 1,1TCOM) 

351 CONTINUE 

C MEASURED SEDIMENT DISCHARGE AT SECTIONS B-l, B-5M, AND B-17 
DO 350 11=1,3 
DO 888 IT=1 , ITCOM 
QLOG=ALOGi 0 (QMES (IT, 1.1 ) ) 

IF(Il.EQ.i) EXPONT=3.05785*QLOG-9 33462 
IF(II .EQ.2) EXPONT=l.92638*QL0G-3.29409 
IF(II.EQ.3) EXPONT = 1.<16664*G)L0G-S 0 3701 
GTOTM(IT)=10.0**EXPONT 
888 CONTINUE 

READ(4,401) NINT 
WRITE(5,401) NINT 

READ(4,402) < INT (I) , G|S< I > , GS (I > , ( PCF (I , J ) , J = 1 ,NFR ) ,1 = 1 ,NINT> 
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WRITE<5,402) <INT<I),QS(I),GS(I),<PCF(I,J),J = 1,NFR),1 = 1,NINT) 
12=11 

DO 309 1=1,NINT 
PCP = 0.0 

DO 301 J = 1 ,NFR 
PC=PCF (I,J)-PCP 
PCFR(I,J)=PC/100.0 
PCP=PCF(I , J) 

301 CONTINUE 
309 CONTINUE 

DO 302 IF=i,NFR 
DO 303 1=1,NINT 
IF(I.GT.l) GO TO 304 
INT1 = INT (1) 

DO 305 J=i,INTi 

305 GMES(J,IF,II)=PCFR (l,IF)*GTOTM(J) 

KC=INT <1) 

GO TO 303 

304 IF(I EQ.NINT) GO TO 306 
IB=KC+i 
IE = INT <I) 

GD=(PCFR(I,IF)-PCFR(1-1,IF))/FLOAT(IE-KC) 

DO 307 J=IB,IE 

307 GMES< J ,IF,I1) = (PCFR (I-1,IF )+GD#FLOAT ( J-KC)) #GTOTM( J ) 

KC=INT <I) 

GO TO 303 

306 CONTINUE 
INT2=INT(NINT) 

DO 308 J=INT2,ITCOM 

308 GMES(J, IF , 11)=PCFR(NINT,IF)*GTOTM(J) 

303 CONTINUE 

302 CONTINUE 
350 CONTINUE 
C 

C BEGINNING OF CHANNEL-BLOCK LOOP 
DO 201 1=1,NSEG 
INL=4 0 

READ(1,403) INLS,SLN 
WRITE(5,403) INLS,SLN 
C GEOMETRIC INPUT FOR CHANNEL BLOCK 

READ(1,404) (XSI(IC),SLOPE<IC),WIDTH(IC>,CPER(IC),EPER(IC), 
&CDEP(IC),EDEP(IC),IC=i,INLS) 

WRITE(5,404) (XSI(IC),SLOPE(IC),WIDTH(IC),CPER(IC),EPER(IC>, 
&CDEP(IC),EDEP(IC),IC=1,INLS > 

DO 801 IC=1,INLS 
PCP=0.0 

DO 802 IF=1,NFR 
PCF(IC,IF)=PCFF(I,IC, IF) 

PCBI(IC,IF)=PCF(IC,IF)-PCP 
PCP=PCF(IC,IF) 

802 CONTINUE 
801 CONTINUE 
C 

C INITIAL AND BOUNDARY CONDITIONS FOR CHANNEL BLOCK 
SLOP=0.0 

DO 11 IT=i,ITCOM 
IF(I-2)3,4,S 

3 QUP(IT)=QMES(IT,i) 

GO TO 6 

4 QUP(IT)=QMES(IT,2)+Q(IT,1) 


J.7A 





GO TO 6 

5 QUP <IT)=Q(IT,2> 

6 QLAT(IT) = 0 . 0 
DO 11 IF=i,NFR 
IF<I-2>7,8,9 

7 CUP <IT,IF)~GMES<IT,IF,1)/<QUP(IT)*GAMA*SPGR> 

GO TO 10 

8 CUP(IT,IF)=(G(IT,1,IF>+GMES1IT,IF,2 >/<GAMA#SPGR)>/QUP <IT) 

GO TO 10 

9 CUP <IT,IF)=G<IT,2,IF)/QUP<IT> 

10 GL.AT (IT,IF)=0.0 

11 CONTINUE 

DO 13 IC = i,INLS 
BEDELV (IC) =BEDLVL (I , IC) 

DO 12 IF=i,NFR 

Cl(IC,IF >=CII(I>IC,IF > 

IF(CI(IC,IF> .EQ.0.0)CI(IC,IF)=CUP(1,IF> 

12 PCBC(IC,IF)=PCBI(IC,IF) 

SLOP=SLOP+SLOPE<IC) 

13 CONTINUE 
SLOP=SLOP/FLOAT<INLS) 

QBASE=QUP(1 ) 

DO 300 IC=i , INL 
QI(IC,1)=QUP<1) 

QI < IC> 2) =Ql)P < 1) 

XI < IC > =Sl.N#FLOAT < IC-1 >/FLOAT (INL) 

300 KSI<IC)---0 

C LIST INITIAL BED ELEVATION AND SIZE COMPOSITION 
URITE(S,50i) I 

URITE<5,502) <IC,<PCF<IC,IF),IF=1,NFR>,BEDELV(IC)>IC=1,INLS) 
C 

C ROUTE UATER AND SEDIMENT THROUGH CHANNEL BLOCK DURING 
C CURRENT TIME STEP 

DO 101 IT = 1,ITCOM 
C 

C UATER ROUTING 
CALL UR0UT2 
C 

C SEDIMENT ROUTING 
CALL SR0UT2 
C 

IF<INL.EQ.0) GO TO 129 

C RESET INITIAL CONDITIONS FOR FLOW CALCULATIONS 
DO 132 J = 1 ,1 Nl. 

XI(J>=XC(J) 

QI(J,1)=QC(J,1) 

QI < J, 2) = QC ( J, 2) 

KSI< J)=KSC < J) 

132 CONTINUE 
12.9 CONTINUE 

C RESET INITIAL CONDITIONS FOR SEDIMENT CALCULATIONS 
DO 270 IF=1,NFR 
DO 271 IC=i , INLS 
IFUC.EQ. 1) GO TO 272 
CI<IC,1F)=CC(IC,IF) 

GO TO 271 

272 CI(IC,JF) =-CUP < IT , IF ) 

271 CONTINUE 
270 CONTINUE 

IF1IT.NE.170) GO TO 101 
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WRITE(5,S0 3) IT 
DO 900 IC=i ,INLS 
PCF(IC,i)=PCBC(IC,l> 

DO 901 IF=2,NFR 

901 PCF(IC,IF)=PCF(IC,IF-1)+PCBC(IC,IF) 

900 CONTINUE 

C LIST BED ELEVATION AND SIZE COMPOSITION AT TIME STEP 170 

WRITE(5,S02> (IC,(PCFtIC,IF),IF=i,NFR>,BEDELV<IC),IC=i,INLS) 
101 CONTINUE 
C DISCRETIZED OUTFLOW 
QP=QBASE 

DO 315 IT=1,ITCOM 
QCC=Q(IT,1) 

Q(IT,I)=(QP+QCC)/2.0 
QP=QCC 

31S CONTINUE 
C 

WRITE < 5 f 504) 

DO 911 IC=1,INLS 
PCF(IC,i)=PCBC(IC,l> 

DO 912 IF = 2,NFR 

912 PCF(IC,IF)=PCF(IC,IF-1>+PCBC(IC,IF) 

911 CONTINUE 

C LIST FINAL BED ELEVATION AND SIZE COMPOSITION 

WRITE(5,5021 (IC, (PCFdC, IF) , IF=1, NFR >, BEDELV (IC), IC=1, INLS > 
C 

DO 913 IC=i , INLS 
BEDLVL (1,10= BEDELV (IC) 

DO 914 IF=1,NFR 
PCFF(I,IC,IF)=PCF(IC,IF) 

CII (I ,IC,IF)=CI(IC,IF) 

914 CONTINUE 

913 CONTINUE 
201 CONTINUE 

C END OF CHANNEL-BLOCK LOOP 
REWIND 1 

C COMPUTE HYDROGRAPH, SEDIMENTGRAPH, AND SEDIMENT YIELD AT 
C SECTION B-17 
TGOUT=0 0 
TGMES=0 0 
DO 299 IF = 1 ,NFR 
TGO(IF)=0 0 
TGM(IF)=0.0 
TGIN<IF)=0.0 
299 CONTINUE 

WRITE(5,505) 

DO 41 IT=1,ITCOM 
QOUT <IT)=Q(IT,NSEG) 

GTOTC(IT)=0.0 
DO 42 IF=i,NFR 

G<IT,NSEG,IF)=G(IT,NSEG,IF)#GAMA#SPGR 
GTOTC(IT)=GTOTC(IT)+G(IT,NSEG,IF> 
TGO(IF)=TGO(IF)+G(IT,NSEG,IF)*DTS 
TGM<IF)=TGM(IF)+GMES(IT,IF,3)*DTS 

TGIN(IF)=TGIN(IF)+(GMES(IT,IF,1)+GMES(IT,IF,2>)*DTS 
42 CONTINUE 

TGOUT=TGOUT+GTOTC(IT)*DTS 
TGMES=TGMES+GTOTM(IT)*DTS 

C COMPUTE AND LIST CHANGES IN SIZE COMPOSITION OF SEDIMENT 
C LOAD AT B-l7 DURING CURRENT TIME BLOCK 





990 


99 i 


992 


43 


902 

41 

C 


999 

C 

C 

333 

401 

402 

403 

404 

405 

406 

407 

408 

409 

501 

502 

503 

504 

505 

506 

509 

510 


511 


IF(N.NE.i) GO TO 990 

IF(IT.EQ.213.OR.IT.EQ.214.OR.IT.EQ.264) GO TO 43 
GO TO 41 

IF(N .NE.2 > GO TO 991 

IF< IT.EQ.3 OR IT.EQ.4 OR.IT EQ.54 OR.IT EQ.99) GO TO 43 
IFdT.EQ.100 OR IT.EQ.144.0R.IT.EQ.148.0R.IT EQ.149) GO TO 43 
IF(IT.EQ.189.OR IT.EQ.191.0R.IT EQ.193.OR.IT.EQ.239) GO TO 43 
IFCIT.EQ 240 OR.IT.EQ.241.OR.IT.EQ.244.OR . IT.EQ.247) GO TO 43 
GO TO 41 

IF<N.NE.3) GO TO 992 

IF(IT EQ 21 OR.IT.EQ 22 OR.IT.EQ.24 OR.IT EQ.69) GO TO 43 
IF<IT.EQ.70.OR IT.EQ.72.OR.IT.EQ.119.OR.IT.EQ.120> GO TO 43 
l!'<IT.EQ.165.OR IT EQ.167.OR.IT.EQ.212.OR.IT.EQ.213) GO TO 43 
IF(IT.EQ.262.OR.IT.EQ.263.OR IT.EQ.264) GO TO 43 


GO TO 41 

IFdT.EQ 57.OR . IT . EQ 100 OR IT .EQ 144 OR . IT EQ 185) GO TO 43 
IF(IT.EQ.186 OR IT.EQ 192 OR.IT EQ.237 OR.IT.EQ.238) GO TO 43 
GO TO 41 
CONTINUE 

PCW(1)=G(IT,NSEG,1)/GTOTC(IT)*100 . 0 
DO 902 IF = 2,NFR 

PCM <IF)=PCW <IF-i)+G<IT,NSEG,IF)/GTOTC <IT)#10 0.0 
WRITE<5,506) N,IT,(PCW(IF),IF = 1,NFR) 

CONTINUE 

LIST SEDIMENT YIELD, HYDROGRAPH, AND SEDIMENTGRAPH AT 8-17 
WRITE < 5,509) TGMES,TGOUT 
WRITE<5,510) 

URITE<5,511) <IF,TGO<IF),TGM(IF),TGIN(IF),IF=i ,NFR ) 

WRITE(5,512) 

WRITE<S,513) <1,(QMESII, J) ,J=i,3), QOUT(I),GTDTM(I),GTOTC(I), 
41=1,ITCOM) 

CONTINUE 

END OF TIME-BLOCK LOOP 


FORMAT(10F7.2) 

FORMAT <14) 

FORMAT<14,F6.1,F8 3,10F6.2) 

FORMAT <14,Fi0.2) 

FORMAT < F7.1,F3.5,F7.3,F7.3,3F6.3) 

FORMAT<10F7.2) 

FORMAT(14,F7.1,14,14) 

FORMAT(3F7.2,FiO.7) 

FORMAT <10F7.2) 

FORMAT < 4F10.5) 

F0RMAT(//15X,'INITIAL SIZE DISTRIBUTION OF BED MATERIAL AND', 

6' BED ELEVATION FOR EACH SECTION OF CHANNEL BLOCK',14/) 

FORMAT <15X,14,10F8.3,F15.7 , ' FEET') 

F0RMAK/15X,'SIZE DISTRIBUTION OF BED MATERIAL AND BED ELEVATION 
4' FOR EACH SECTION OF CHANNEL BLOCK AT TIME STEP',15/) 
F0RMAT(/15X,'FINAL SIZE DISTRIBUTION OF BED MATERIAL AND BED', 

4' ELEVATION FOR EACH SECTION OF CHANNEL BLOCK'/) 

F0RMAT<//15X,'VARIATION OF LOAD SIZE DISTRIBUTION AT B-17 ', 

4' DURING CURRENT TIME. BLOCK'/) 

F0RMAT(5X,'TIME BLOCK',12,3X>'TIME STEP',14,SX,10F8.3) 
F0RMAT<//15X,'TOTAL MEASURED SEDIMENT YIELD =',F20.S,' LBS.'/ISX 
4'TOTAL COMPUTED SEDIMENT YIELD =',F20 5,' LBS.'//) 

F0RMAT(//13X,'SEDIMENT COMPUTED MEASURED 

4,' MEASURED'/13X,'TRACTION YIELD(LBS) YIELD', 

4'(LBS) INFL OW(LBS)'/) 

FORMAT(15X,13,F20.S,F20.5,F20 5> 


J. 7 7 





512 F0RMAT<//S4X,'MEAS.'/16X,'TIME COMPUTED FLOW(CFS) AT 

&, 'FLOW SED. LOAD(L.BS) AT B-17'/16X, 'STEP B-l B-SM', 

t.' B-17 B-17 MEAS. COMP.'//) 

513 FORMAT <15X,14,6F10.3) 

520 FORMAT(2X,F20.5/2X,E15.7) 

END FILE 5 
STOP 
END 
C 

SUBROUTINE WR0UT2 
C 

C THIS SUBROUTINE ROUTES WATER THROUGH A CHANNEL BLOCK USING THE 
C KINEMATIC WAVE SCHEME DEVELOPED BY BORAH ET AL<1980). A LIST OF THE 
C IMPORTANT VARIABLES USED IN THIS SUBROUTINE IS GIVEN BELOW: 

C - 

C NAME DESCRIPTION UNITS 

C - 

C IC 

C 

C INL 
C 

C IS 

C 
C 

C KSC(IC) 

C 
C 
C 

C KSI(IC) 

C 
C 
C 

C Q(IT,I) 

C 
C 

C QC(IC,L) 

C 
C 

C QI (IC , L.) 

C 
C 

C XC(IC) 

C 
C 

C XI(IC) 

c 
c 

c - 

c 

COMMON /ROUT/ I , IT ,SLN ,CHEZY ,DTS , ITC.OM , INL ,QBASE , 

&QUP(300),QI(50,2),XI<SO ),KSI(50),Q<300,3) 

COMMON /WROUT/ SLOP,QLAT<300>,QC<50,2>,XC(50>,KSC<50> 

EXP-1 S 
BET=1.0/EXP 
EXP 1=EXP +1.0 
EXMl=EXP-i.0 
KIN=CHEZY*SLOP**0.5 
TF.RM=EXP*KIN*DTS 
QL-QLAT <IT) 

QU=QUP(IT) 


PARAMETER USED TO LABEL INDIVIDUAL WAVES. 

NUMBER OF NODE POINTS IN THE CHANNEL BLOCK 

COUNTS NUMBER OF SHOCK WAVES OCCURRING AT THE END 
OF THE CURRENT TIME STEP . 

FLAG USED TO CHARACTERIZE THE WAVE IC OCCURRING IN 
THE CHANNEL BLOCK, AT THE END OF THE CURRENT TIME 
STEP. KSC=0 FOR CHARACTERISTICS, KSC>0 FOR SHOCKS. 

FLAG USED TO CHARACTERIZE THE WAVE IC OCCURRING IN 
THE CHANNEL BLOCK, AT THE START OF THE CURRENT TIME 
STEP. KSI=0 FOR CHARACTERISTICS, KSI>0 FOR SHOCKS. 

COMPUTED WATER OUTFLOW FROM THE CHANNEL BLOCK I , AT C.FS 
THE END OF THE CURRENT TIME STEP IT. 

FLOW DISCHARGE AHEAD <L=i> OR BEHIND <L=2) OF THE CFS 
SHOCK IC, AT THE END OF THE CURRENT TIME STEP. 

FLOW DISCHARGE AHEAD (L=l) OR BEHIND <L=2> OF THE CFS 
SHOCK IC, AT THE START OF THE CURRENT TIME STEP. 

DISTANCE OF WAVE FRONT IC TO UPSTREAM END OF THE FT 

CHANNEL BLOCK, AT THE END OF THE CURRENT TIME STEP. 

DISTANCE. OF WAVE FRONT IC TO UPSTREAM END OF THE FT 

wHANNEL BLOCK, AT THE START OF THE CURRENT TIME 
STEP . 








IMQL.EQ. 0 . 0 ANP.RU.FR. 0 . 0) GO TO 130 
C PROJECT ALL CHARACTERISTICS TO NEW TIME LEVEL 
AC=(QU/KIN)**BET+QL*DTS/2 0 
QC(i,l)=KIN*AC**EXP 
IF < QL.EQ 0.0) GO TO 102 
XC<l)=<QC<i,i)-QU>/QL 
GO TO 103 

102 XC<i>=TERM*AC**EXMi/2 0 

103 ICST = 1 
KSC(1>=0 
GO TO 131 

130 ICST=0 

131 IF<INL.EQ.0) GO TO ISO 
IS=0 

DO 104 IC=1,INL 
IA=IC+ICST-IS 

IF(KSKIC) -EQ.0 OR QI<IC,1) .GE.QI(IC,2>> GO TO 10S 
C PROPAGATION OF SHOCK WAVE 
AA=(QI(lC,i)/KIN)**BET 
AB-< QI<IC,2)/KIN)**BET 
AAF=AA+QL*PTS 
ABF = AB+QL.*DTS 
QC< IA,1 )=KIN#AAF**EXP 
QC<IA,2)=KIN*ABF**EXP 
IF(QL.EQ.0.0) GO TO 10B 
PROD=ALP/<EXPl*<AB-AA)#QL> 

XC<IA)=XI<IC>+PROD*<ABF**EXPl-AAF*#EXPl-AB**EXPi+AA**EXPi) 
GO TO 107 

108 XC<IA)=XI<IC) + <QI(IC,2)-G)I(IC,1> >*DTS/(AB-AA> 

GO TO 107 

C PROPAGATION OF CHARACTERISTIC WAVE 

105 KSI<IC)=0 

AC=(QI (IC , 1>/KTN>**BET+QL*DTS 
QC(IA,1)=KIN*AC*#EXP 
IF < QL.EQ.0.0) GO TO 106 
XC<IA)=XI(IC) + (QC(IA,1) - QI(IC,1))/QL 
GO TO 107 

106 XC<IA)=XI<IC)+TERM*AC**EXM1 
C CHECK FOR NEW SHOCK FORMATION 

107 IF<KSI<IC)GT.0) GO TO 10? 

IF<IA.EQ.1> GO TO 140 
IF<XC(IA).LE.XC(IA-l)) GO TO 110 

C NO SHOCK IS FORMED 
140 KSC(IA)=0 
GO TO 104 

C SHOCK IS FORMED 
110 IA=IA-l 

IF < K SC(IA).GT.0) CO TO 112 
C SHOCK IS FORMED BY TWO CHARACTERISTIC WAVES 
XC< IA) = (XC.(1A)+XC( IA+i ) )/2.0 
QC<IA,1>=QC<IA+i,1) 

QC(IA,2)=QC(IA,1) 

IS--=IS+i 
K SC < IA ) = 1 
GO TO 104 

C THE CHARACTERISTIC WAVE JOINS THE SHOCK AHEAD OE THE FRONT 
112 XC < IA ) =XC ( IA) 

QC(IA,1)=QC<IA +1,1) 

IS=IS+i 
K SC<IA)=11 
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GO TO 104 

C CHECK IF THE PROPAGATING SHOCK IS INTERSECTED BY ANY OTHER SHOCK 
C OR CHARACTERISTIC WAVE 
109 IF <IA EQ.1) GO TO 141 

IF(XC<IA> LE.XC(IA-D) CO TO 113 
C THE PROPAGATING SHOCK IS NOT INTERSECTED 
141 KSC < IA ) = 1 

GO TO 104 

C THE PROPAGATING SHOCK IS INTERSECTED 

113 IA=IA-i 

IF <KSC<IA>.GT.0) GO TO 114 

C THE PROPAGATING SHOCK IS JOINED BY A CHARACTERISTIC WAVE 
C BEHIND THE FRONT 

XC<IA)=XC(IA+1) 

OC<IA,2)=QC<IA , 1 ) 

QC<IA,1)~QC<IA+1,1) 

ISdS+1 
KSC(IA > ~11 
GO TO 104 

C NEW SHOCK IS FORMED BY TWO INTERSECTING SHOCKS 

114 XC(IA)=<XC<IAI+XCCIA+1)>/2 0 
QC<IA,i)=QC<IA+i ,1) 

QC(IA,2)=QC(IA,2) 

IS=IS+1 
KSC <IA)=2 
104 CONTINUE 

C COMPUTE OUTFLOW FROM CHANNEL BLOCK DURING CURRENT TIME STEP 
ISO CONTINUE 
IAD=0 

IF(INL.EQ.O) IA=ICST 
IF< IA . Et). 0 ) GO TO 122 
XB=0.0 
QB=QIJ 

DO 115 J = 1 , IA 

IF(XC(IA) .GE.SI.N) GO TO 117 
XB---XC(IA) 

QB=QC(IA,1) 

IF<KSC(IA). GT . 0 ) QB=(QC(IA,l)+QC<IA,2>>/2.0 

117 IC=IA-J+1 

IF(XC<IC) LT.SLN) GO TO 115 

IAD=IAD+1 

XA=-XC(IC) 

QA-QC(10,1) 

IF < KSC <IC) .GT.0) QA= < QC(IC,1)+QC<IC,2>)/2.0 

115 CONTINUE 

IF(IAD.GT 0) GO TO 118 
AIMQBASE/KIN)**BET 

IF<IT.GT.1) AI=(Q(IT-1,I)/KIN)**BET 

AC=AI+QL*DTS 

QA=KIN*AC**EXP 

XA=TERM*AC*#EXM1+SLN 

IF <QL GT . 0 . 0) X A= < QA-Q (IT -1 , I ) )/Qt„+SLN 

118 Q(IT,I)~QB+< QA-QB)* < SLN-XB >/(XA-XB > 

INL=IA-IAD 

GO TO 123 

122 Q(IT >I ) = 0 . 0 
INL = 0 

123 RETURN 
END 

C 
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SUBROUTINE SR0UT2 


THIS SUBROUTINE PERFORMS SEDIMENT CALCULATIONS AT ALL NODE POINTS 
IN THE CHANNEL BLOCK DURING CURRENT TIME STEP 

DIMENSION WATMAT(10 >,BEDM AT(10 ) ,SCAP(10>,ERS(10) ,ACCh<i0,10 ) , 

4>TCA( 10 > 

COMMON /ROUT/ I,IT,SIN,CHEZY,DTS,ITCOM,INL,QBASE, 

&QUP(300),QI(50,2),XI(SO),KSI(SO),Q(300,3) 

COMMON /SROUT/ SPGR,GAMA,SNU,NER,INLS,WIDTH< i 0),SLOPE(10 ) , 
iCPER<10),EPER <10 > ,CDEP < 10 ) ,EDEP ( 10 > ,GLAT < 3 00,10 ) ,XSI ( 10 ) , 
iXSCCl0,10),BEDELV(10),CI(i0 ( 10),CC(10,10),PCBI<10,10),PCBC(10,10), 
&CUP(300,10),DMM(10),COEF(10),OS(10),POR(10),G(264,3,10),PCU(10 >, 
4ER0 ,CHI,CEL 

DO 204 IC=i , INI.S 

INTERPOLATE FLOW VALUES AT THE NODE POINTS 
IFdC.EQ.l) GO TO 20S 
IF(XSI(IC> L.T XI(1)) GO TO 206 
DO 207 J = 1 ,100 

IF(XSI(IC).GE,XI(J) AND.XSI(IC).LT.XI(J+i)) GO TO 208 
GO TO 207 

208 Ql=QT(J , 1) 

IF(KSI(J) GT.0) Ql=(QI(J,i)+QI(J,2))/2.0 
Q2 = QI <H 1,1) 

IF < KSI(J +1) . GT 0 ) Q2=(QI(J + 1,1)+QI(J+2,2))/2.0 
QP = Ql + (Q2-Qi)*<XSIdC)-XI( J))/(XI(J + l)-XI(J) ) 

GO TO 209 
207 CONTINUE 
GO TO 209 

206 IFdT.EQ. 1) GO TO 210 
Q1=QUP(IT-1 ) 

GO TO 211 

210 Qi=QI(IC,1) 

211 Q2=QI(l,i) 

IF(KSI(1).GT.0 > 02-(QI(1,1)+QI(1,2)>/2.0 
RP-QI+ < Q2-Q1)*XSI<IC)/XI(1) 

GO TO 209 

20S IFdT.EQ.1) GO TO 212 
G!P=G)UP (IT-1 ) 

GO TO 209 

212 QP = RI (IC , 1 ) 

209 CONTINUE 
C 

C UPDATE THE HYDRAULIC PARAMETERS AT THE NODE POINTS 
SLP=SL OPE(1C) 

UH)TH=UIDTH( IC) 

CPR-CPER (JO 
F-PR-EPER(IC) 

CDP =CDEP(IC) 

EDP-EDEP(IC) 

E X P = 1 . S 
BE T = 1 0/EXP 
KIN=CHF.ZY*SL P**0,5 
QE-QP 

AE- ( QE/K IN ) **EiET 

IF (AE . L.T 1 0E-5) GO TO 213 

VLL-QE/AF 

PPTH=CDP*AE**F..DP 
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WF.P=CPR*AE*#EPR 

RHR=AE/UEP 

HYR=DPTH 

UST“SQRT(32.2*HYR*SLP ) 

IF(IC.EQ.INLS) GO TO 241 
DELX=XSI(IC+1)-XSI(IC> 

GO TO 242 

241 DELX=SLN-XSI<INLS) 

242 CONTINUE 
C 

C COMPUTE TRANSPORT CAPACITIES AND TERM DELTA (EQ. 11) 

SUM = 0.0 

DO 220 I F = 1,NFR 
GL=GLAT(IT,IF) 

CP=CI<IC,IF) 

D = DMM(IF) 

DFT=DMM(IF)/304.8 
U=VS(IF) 

IF<DMM(IF),LE.4 0) GO TO 860 
IF ( DMM( IF ) . l.T . S . 0 )GO TO 861 
MEYERP=8.0 

CALL MEYER(MEYERP,UST,RHB,DFT,GAMA,SNU,SPGR,WDTH,QE,VEL,CONC) 
SCAP <IF)=CONC*AE/l.0E6/SPGR 
IF < DhM(IF) GT.30.0) SCAP<1F)=0.0 
GO TO 862 

860 CALL YANG(DFT,UST,SNU, VEL,SLR,U,SPGR, CONC) 

SCAP(IF)=CCjNC*AE/1.0E6/SPCR 

GO TO 862 

861 CALL DUBOYS(D,DFT, AE,WDTH,SLP,SPGR,SNU,GAMA,CHI,UST,SFLOU> 
SCAP(IF)-SFLOW/VEL 

862 WATMAT < IF ) ==CP*AE + GL*I)ELX/VEL 
IF(SCAP(IF) EQ.0.0) GO TO 220 
5UM=SUM+WATMAT(IF)/SCAP(IF) 

220 CONTINUE 
C 

C COMPUTE ACTIVE (OR ARMOR) LAYER THICKNESS 
PARM= 0 0 
NARM=NFR 
DARM=DMM(NFR) 

DO 810 IF =-1 , NFR 
IF(SCAP(IF> GT.0.0) GO TO 810 
IF(IF EQ.1) GO TO 811 
IF(SCAPdF-i) .GT.0.0) CO TO 811 
GO TO 812 

811 NARM~IF 
DARM=DMM(IF) 

812 PARM=PARM + PC.E<C( IC, IF ) 

810 CONTINUE 

IF(FARM.CT 0.0) GO TO 814 

DO 81S IFF = 1,NFR 

IF A = NFR-IFF f1 

NARM=IFA 

DARM=DMM(IFA) 

PARM = PCK<C(IC,IFA) 

IF(PARM.CT.0.0> GO TO 814 
81S CONTINUE 

814 ARMHT =10 0 (UDARM/PARM/304 8 

C COMPUTE VUL1JMF OF RED MATERIAL FRACTIONS IN ACTIVE LAYER 
DO 813 IF=1,NFR 

REDMA 1 (IT ) ARMHT*PCRC(IC,IF>*WDTH/100.0 


813 





IF<1.0-SUM) 222,223,224 

continue; 

DEPOSITION LOOP 
TBM-0 . 0 
BEDUP=0 0 
DO 230 J=1, NFR 
IF = NFR-J +1 
CP=CI<IC,IF) 

IF(SCAP(IF) . EQ . 0 0) CO TO 231 
RESCAP=SCAP<IF)*<1.0-SUM) 

DF.PO=0 0 

IF CRESCAP GE 0 0 ) GO TO 232 
IFCABSCRESCAP ) . GT UATMAI ( IF) ) C,0 TO 231 
DEPO=ABS(RESCAP) 

GO TO 232 
DEPO-WATMAT(IF ) 

BT-2 0*VS<1F)*DTS/RHB 
IF(BT.L T 1 0)DFPO-BT*DEPO 
UJAT MAT ( IF ) - UATMAT < IF)-DEPP 
ADI) = 0 . 0 

IF(SCAP(IF) GT.0 0) ADD- DEPO/SCAP(IF) 
SUM-SUM-ADD 

UPDATE LOAD AT NODE POINTS 
CAP-SCAP <IF)/AE 
CCC=UATMAT<IF>/AE 
XKP=CAP-0 999*CP 
XK C = CAP-0 999*CXC 
ADD-0 . 0 

IF (CAP GT 0.0) ADD'CLL *f.AP / ( AE*XKP *XKC ) 

IF<ADD LT 0.0) ADD-0 0 
FNLR=1 0+ADD 
DELT=DELX*FNLR/VEL 
IF < DELT.LT.DTS > GO TO 220U 
IFUC.EQ. INLS) GO TO 233 
CCI-CI<IC+1,IF) 

CC<IC+1,IF) = CCI+<CCC-CC1 )*i>TS/BFL I 
GO TO 234 

SEDIMENT OUT FLOW FROM CHANNEL BLOCE: 

IF(IT EQ 1) CO TO 23G 
CCI-G<IT-1,1,IF)/0<IT-1 ,1 ) 

GO TO 236 
CCI-CI<IC,IF) 

G < I T , I , IF ) -' (F.GI + ( CCC-CCI ) *I)TS/I>f L T ) #0 < IT , I > 

GO TO 234 

0 XDEL-DLLT*VKI../I N! R 

IF <IC EQ.1) GO TO 2201 
CC,l! = CC(IC, IF) 

GO TO 2202 

1 CCB~CUP(IT,IF) 

2 CC(IC+1,IF)-CCB+(CCC-CCEO*DELX/XDEL 
IFUC.F.Q INLS > G< IT , I , IF >-CCt JL + ) , IF >*Q< IT , I > 
BFDMAT(IF)-BEDMAT(IF > +DEPO*DTS/DrLT 

TBM = TBM + BE DMA1<IF ) 

BEDUP -BF.DUP + DEPO'#DTS /DELT/UDTH/POR < IF ) 

CONTINUE 

CHANGE IN BED ELEVATION AND BED LAYER COMPOSITION 
BEDELV( IC ) =BE DI.LV < IC)+BE DUP 
DO 239 IF=1,NF R 

PCBC <IC,IF)=BEDMAT <IF)/EPM*100 . 0 










GO TO 204 

223 CONTINUE 
C 

C EQUILIBRIUM LOOP 

DO 240 IF= 1,NFR 
CP=CI (IC , IF') 

CCC=WATNAT(IF)/AE 
CAP =SCAP(IF>/AE 
XKP=CAP-0 999*CP 
XKC=CAP-0 9?9*CCC 
ADD=0 0 

IF(CAP .GT.0 . 0) ADD=CEL*CAP/(AE#XKP#XKC> 

IF<ADD.LT 0 0) ADD-0 0 
FNl.R = 1 0+AI)D 
DELT+DELX*FNLR/VEL 
C UPDATE lOAD AT NODE POINTS 

IF(DELT L T DTS) CO TO 2300 
IF< IC EQ INLS) GO TU 243 
CCI=C1(IC+1,IF) 

CC( IC+1,IF >~CCI + (CCC-CCI)*DTS/DCLT 
GO TO 240 

C SEDIMENT OUTFLOW FROM CHANNEL BLOCK 
243 IF(IT EQ 1) GO TO 244 

CCI=G(IT-i,I,IF)/Q(IT-i , I ) 

GO TO 24S 

2*4 CCI-CI(IC,IF) 

24S G(IT,I,IF)=(CC1+(CCC-CCI)*DTS/DELT)*QCIT,1) 

GO TO 240 

2300 XDEL-DLLT fcUEL/FNLR 

IF(IC EQ.1) GO TO 2301 
CCB=CC(IC,IF) 

GO TO 2302 

2301 CCB-CUP(IT,IF) 

2302 CC<IC+1,IF)=CCB+(CCC-CCB>*DELX/XPEL 
IFUC.FQ I NLS) G(IT,I,IF)=CC(IC + 1,IF)#Q(IT,I) 

240 CONTINUE 
GO TO 204 

224 CONTINUE 
C 

C EROSION AND ARMORING CALCULATIONS 
C 

C COMPUTE VOLUME ENTRAINMENT MATRIX 
DO 201 I FA-1 , NFR 
DEN=0 0 

LOOP-NFR-IFA+1 
DO 282 J = 1.LOOP 
IF-J+IFA-1 

282 DEN-DFN+COEF(J)*BEDMAT(IF ) 

DO 283 J=1.LOOP 
IF-J+IFA-1 

IF(DEN EQ 00) GO TO 030 

ACCM(IF,TFA)=COEF(J)*BEDMAT<IF >/DEN*BEDMAT(IFA) 
GO TO 283 

830 ACCM<IF,IFA)=0.0 

283 CONTINUE 
281 CONTINUE 
C 

C EROSION LOOP 

DO 270 IF =1,NFR 
TCA(IF)=0.0 
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DO 279 IFA=i,IF 

279 TCA (IF)-TCA(IF)+ACCM(IF,IFA> 

278 CONTINUE 

DO 237 I F =1,MFR 
237 ERS<IF)=0.0 

DO 284 IF-1,NFR 

IF < SCAP(IF) , EQ 0.0) GO TO 284 

RESCAP=SCAP(IF>*(1 0-SUH) 

IF < RESCAP LE O 0> GO TO 285 

IF C RESCdP .GE.TCA(IF)) GO TO 288 

DO 831 J = 1 , IF 

DO 287 Jl=--i , IF 

IF < J1.EQ.1) IF A =IF 

IF(Ji.GT.l) IFA-Ji-i 

IF < SC AP < I F A ) FQ 0 0) GO TO 287 

resc:ap=scap ( ifa)*( i o-sum> 

CNTBN=ACCM ( IF , IF A)/FLOAT ( IF > 

IF < RESCAP GT.CNIDN) GO TO 280 
ER OSN-ERO*R I.SCAP 
GO TO 289 

288 EROSN^I RO*CNTBN 

289 ERG ( IFA)'ERS( II A)+EROF>N 
Sim=-SUM+EROSN/SCAP ( IFA) 

IF(SUM CF 1.0) GO 10 288 

287 CONTINUE 
831 CONTINUE 

288 CONTINUE 

DO 290 IFA- 1 ,1F 

IF(SCAP(IF A) IQ 0 0) CO TO 290 
EROSN ERO*ACCM<IF, IFA > 

ERS<IFA)=EKS(II A)+FROON 
SUM-SUMtER0SN/5CAP(IFA) 

290 CONIINUI 

284 CONTIMUt 

C UPDATE LOAD A I NODI POINTS 

285 CONTINUE 
TBM - 0 0 
DEDUN-O 0 

DO 2S0 IF-).NFR 
CP Cl ( IC. , I F ) 

UA1 MAI (IF; WA1MAT ( IF ) t[.RS( IF > 

CCC=UA1 MAT ( I F ) /,*,£ 

CAP = SC AH ()t )/AE 
XKP - CAP - 0 999♦CP 
XKC-CAP 0 999HCCC 
ADD-0 0 

IF" ( CAP GT . 0 0) AI>P~CF I *CAF>/(AF.*/KP*XKC) 

IF (ADD L. 1 0 0 ) ADD- 0 0 

FNL.R = 1 0tADD 

DELT-DElX#F NLR/UEL 

IF (DELT.L.T DTS) GO TO 24 0 0 

IF( IC EQ INLS) GO 10 2S1 

CCI =CI(IG + 1 , IF) 

GC ( IC + 1 , IF ) -CCI + < CCC—CCI >*DT3/I)FLT 
GO TO 252 

C SEDIMENT OUTFLOW FROM CHANNEL BLOCK 
251 IF(IT EQ.l) CO TO 253 

CCI~G(IT-1,I,IF)/Q(IT-1,I) 

GO TO 254 

253 CCI = CI <IC,IF) 
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254 GUT, I ,IF)=(CCI+<CCC-CCI>*DTS/DELT)*Q(IT,1> 

GO TO 252 

2400 XDEL=DELT*VEL/FNLR 
IFUC . EQ.1> GO TO 2401 
CCB=CC(IC,IF) 

GO TO 2402 

2401 CCB-CUP(IT,IF) 

2402 CCUC + l ,IF)=CCB+(CCC-CCB)*DELX/XDEL 
IFUC.EQ. INLS) G ( IT , I, IF) =CC < IC+1, IF > *Q< IT , I > 

252 TERS-ERS(IF)*DTS/DELT 

IF(TERO.GT.BEDMAT(IF)) TERO=BEDMAT(IF) 

BEDMAT (IF)“BEDMAT(IF)-TERO 
TBM=TPM+BEDMAT(IF) 

C UPDATE BED ELEVATION 

BEDWN-BEDWN+TERO/WDTH/PGR(IF) 

250 CONTINUF 

BEDEL V( IC ) BLDFI.V ( IC ) -BEDWN 
C ADJUSTMENT OF ACTIVE l.AYLR THICKNESS 
IF(TBM GT 0 0) CO TO 294 
DO 293 IF - 1 , NFR 

293 PCBC <IC , 1F)=PCBI(IC,IF > 

GO TO 204 

294 CONT1NUF 
PARM-0 0 

DO 295 IF-5 , NFR 

PCBC<IC,IF)=BEDMAT <IF)/TBM*100 0 
IF < SC.AP (IF) GT 0 0) GO TO 295 
PARM = PARM-t PCBC( IC , IF ) 

295 CONTINUE 

IF(PARM Cl 0 0) GO TO 820 

DO 821 1 FF =1 .NFR 

IFA=NFR-IFF+1 

NARM = II A 

DAR M = DMM(If A) 

PARM-PCBC(IC.IF A) 

IF(PARM GT 0 0) GO TO 820 
821 CONTINUE 

820 CURARM*1 0 0 0 /? ARM*DARM/30 4 8 

ADHT-CURARM-TBM/WDT H 
IE(ADHT IE 0 0) GO TO 204 
ADARM-ADHT 
TBM -0 0 

DO 296 IF - 1 .NFR 

BE.DMAT ( IE ) - BFDMAT ( IE ) t ADARM*PCB I ( IC ,IF)/1 00 0*UDTH 

296 TBM=TBM«BEDMAT(IF) 

DO 297 IF -1, NFR 

297 PCBC(IC,IF)-BEDMAT<IF)/TBM*10 0 . 0 
GO TO 204 

213 CONTINUE 

C DEPOSITION CAUSED BY TERMINATION OF FLOW 



DO 260 IF -1 ,NFR 
IFUC.EQ INLS) GO TO 

261 

261 

CCUC+l, IF) = 0.0 

GO TO 262 

GUT,I, IF) =-0 . 0 


262 

BEDMAT(IF)“BEDMAT(IF 

) + CIUC,IF>*AE 

260 

CONTINUE 


204 

CONTINUE 



RETURN 

END 
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S UBROUTINE DUBO YS (D , DFT , AE , WDTH , SLF , SPGR , SNU , GAMA , CHI, 

&UST,SFLOUI) 

THIS SUBROUTINE COMPUTES THE TRANSPORT RATE OF NONCOHESIVE 
SEDIMENTS USING A DUBOYS TYPE FORMULA(GRAF, 1971) 

CALL SHIELD<DFT,SPCR,SNU,GAMA,UST,TC> 

CHI=CHI/(D**0 7S) 

TAO=GAMA*AE/WDTH*SlP 
SFLOU=WDTH*CHI*TAO*(TAO-TC) 

IF < SFLOU.LE.0.0)SFLOW=0.0 

RETURN 

ENI) 

SUBROUTINE MEYER<MEYERP,USTAR,RHP,D,WEIGHT,V T SC,S,UIDTH,Q,V,GB,XO) 

THIS SUBROUTINE COMPUTES THE TRANSPORT RATE OF NONCOHESIVE 
SEDIMENTS USING THE BEDLOAD FORMULA DEVELOPED BY MEYER-PETER 
AND MULLER(1948) 

FCTN<F,RO,REY)-(i /SQRT(F) )-2 . * At OG10 < RO > -1 1-1 + 2 . *ALOG10 < 1 +9 35* 
*RO/<REY*SQRT<F) ) ) 

G = 32.17 

REY-4 0*RHB*V/VISC 

R0=2.*RHB/D 

RSTAR = 2 #D*USTAR/VISC 

IF < RSTAR GT 70.0) GO TO 7 


F 1 = 0 OOfe 
F2=0 09 
1 = 0 

FCT1=FCTN(F1,RO,REY ) 

F3=0 S*<F1+F2> 

FCT3=FCTN(F'3,R0,REY.) 

IF(FCTl*r CT3) 1 ,5,2 
F2-+F3 
GO TO 3 
F1 =F 3 

IF(APS(F 2 -F 1 ) - 0 0001)5,5,4 
1 = 1+1 

IF < I L T 100) GO TO F. 

GO TO 7 
F -F 3 
GO TO O 

F = < 1 0/(2 * Al OC1 0 < R ()) + 1 14))**?. 

AK=V*SQRT(F/H )/USTAP 
AK=0 75 

Y = (USTAR*USTAR )/( (S-l 0 )*C.*D) 

Gi=MEYERP*WIDTH*(UST AR**3 0)*(WEIGHT/G)*(S/(S-1.0)> 
CALL SHIELD CD ,S , VISE. ,UCIGHT ,USTAR ,TC> 

C TSTAR=TC/<<S-1.0)*WEIGHT*D> 

TST AR = 0.047 
G2=AK**1.5-(T STAR/Y) 

IF < G2 . LT . 0 . 0 ) C.0 TO 9 
G?=G2**1.5 
GB=G1*G2 

XO= < GB*1.0E6)/<«*UEICHT) 

GO TO 10 
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9 


10 

C 


GB = 0 0 
XO=0 0 
RETURN 
END 


SUBROUTINE YANGTDFT,UST,SNU,V,SLP,W,S,CONC> 


THIS SUBROUTINE COMPUTES THE TRANSPORT RATE OF NONCOHESIVE SEDIMENTS 
USING THE TOTAL LOAD FORMULA DEVELOPED BY YANG(1973) 


D = DFT 

A=UST*D/SNU 

IF(A GL 70.0) GO TO 7 

VCW=2 5/< ALOGi0 < A>-0.06 > + 0.66 

GO TO 8 

VCW=2.OS 

ESP=V*SLP/W~VCU*SLP 
IFIESP) 9,9,10 
CONC=0 0 
GO TO 11 

0 Fi=S.435-0 286*ALOGi0<U*D/SNU)-0 4S74AL0G10(UST/U) 

FP.- l . 799-0 409*ALOG1 0 <W*D/SNU>-0 314*ALDC10(UST/W) 

F3~AL0G10(ESP) 

E=F1+F2*F3 

c=io o**r 

CONC=C 

1 CONTINUE 

RETURN 
END 

SUBROUTINE SHIELD<D,S,VISC,WIICHT,USTAR,TC) 

THIS SUBROUTINE COMPUTES THE CRITICAL BED SHEAR STRESS DERIVED 
FROM SHIELDS' FUNCTION 
REY - UST AK*l)/VISC 
IF < RFY CT 10 0) CO TO 1 
TC--=0 0S*(S-1 . 0>*WE1CHT*D/REY**0 4 
GO TO 3 

IFIRtY GT 500 (1 ) GO TO 2 

TC*0 022*(S-1 0)*WEIGHT*D*RFY**0 16 

GO TO 3 

TC---0 0 6 * ( S -1 0)*UFICHT*D 

CONTINUE 

REIURN 

END 

SUBROUTINE SETUF.I. ( D ,U ) 

THIS SUBROUTINE COMPUTES THE SETTl 1NG VEIOCITY OF SEDIMENT PARTICLES 
DIMENSION A(2,11) 

DATA A(1,1),A<J,?),A(l,3),A(l,4),A(l,5),A(l,6>,A<l,?>,A<i,8>, 

4A<1,9),A<1,10),A(1,11)/ 

40.04,0 06,0.10,0 J0,0 40,0 80,1 50,200,300,700,10. 0(1/ 

DATA A(2,1),A(2,2>,A<2,3>,A<2,4>,A<2,5>,A<2,6>,A(2,7>,A<2,8), 

4A < 2,9),A < 2,10),A ( 2,11>/ 

40.14,0 32,0.76,2 20 ,5.30,10 50,16.90,20.30,25.60,39.50,44.00/ 

IF<D LE 0 04) GO TO 20 
IF < D GE.10.0) GO TO 21 
GO TO 22 

20 W-0.14/0.04*D 
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21 


GO TO 100 

W=4.5/3.0#(D-i0.0)+44.0 
GO TO 100 
22 CONTINUE 

DO 10 1=1,11 

IF < A<1,1) GT.D) GO TO 11 
Di = A<1,1) 

U1=A<2,I) 

GO TO 10 

11 D2=A<1,I> 

U2=A(2,1) 

GO TO 12 

10 CONTINUE 

12 W=<U2-Uli)*<D-Di)/<D2-Di)+Wi 

100 RETURN 

END 
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